We will start by loading some useful packages. Using the
library function, we can load these packages if they have
previously been installed in your R. If they have not been
previously installed, replace the library function with
install.packages function and place the package name in
quotation marks. That will install the package. Next, run the
library function again and it will load the package. You
only need to install a package once.
library(ggplot2) # Necessary for `ggplot` function.
library(dplyr) # Necessary for `%>%`.
library(kableExtra) # Necessary for `kbl` function.
library(MASS) # Necessary for 'glm.nb' function.
library(car) # Necessary for 'Anova' function.
library(effects) # Necessary for 'effects' function.
library(Rmisc) # Necessary for the 'summarySE' function.
library(pscl) # Necessary for `zeroinfl` function
library(lmtest) # Necessary for 'lrtest' function.
library(betareg) # Necessary for 'betareg' function.
library(emmeans) # Necessary for emmeans function.
library(reshape2) # Necessary for 'melt' function.
library(report) # Necessary for 'report' function
library(ggbeeswarm) # Necessary for geom_beeswarm.
library(ggpubr) # Necessary for ggarrange
library(survival) # Necessary for survfit
library(ggfortify)
library(glmmTMB) # Necessary for 'glmmTMB' function.
library(lme4) # Necessary for mixed lmer
citation("lme4")
## To cite lme4 in publications use:
##
## Douglas Bates, Martin Maechler, Ben Bolker, Steve Walker (2015).
## Fitting Linear Mixed-Effects Models Using lme4. Journal of
## Statistical Software, 67(1), 1-48. doi:10.18637/jss.v067.i01.
##
## A BibTeX entry for LaTeX users is
##
## @Article{,
## title = {Fitting Linear Mixed-Effects Models Using {lme4}},
## author = {Douglas Bates and Martin M{\"a}chler and Ben Bolker and Steve Walker},
## journal = {Journal of Statistical Software},
## year = {2015},
## volume = {67},
## number = {1},
## pages = {1--48},
## doi = {10.18637/jss.v067.i01},
## }
We will set the working directory using the setwd
function. If you are looking to run this code in your computer, be sure
to change the path to your data! Normally we can load dataframes using
the read.csv function.
setwd("~/Desktop/R") # Set the working directory
Fern <- read.csv("Data/FernInterference2023.csv")
DryMass <- read.csv("Data/InterFerDryMass.csv")
WetMass <- read.csv("Data/InterFerWetMass.csv")
Prior to conducting any analyses, we want to make sure that our
variables have the correct class. Often non-numerical variables imported
into R are treated as character variables. We can convert
them into factors using the as.factor function. We can
convert multiple variables into factors by combining
as.factor with the mutate function from the
dplyr package as well as the across function
on a concatenated list of our variables. We will do the
same using the as.Date function to give our date variables
the proper class.
Fern <- Fern %>%
dplyr::mutate(across(c("ID", "Plot", "Species", "Aboveground", "Belowground",
"Germinated", "State"), as.factor))
Fern <- Fern %>%
dplyr::mutate(across(c("Germination.Date", "Death.Date", "Coty.Date",
"Date.First.Leaf", "Date.Highest.Leaf"),
~as.Date(., format = "%d-%b-%y")))
WetMass <- WetMass %>%
dplyr::mutate(across(c("ID", "Plot", "Species", "Aboveground", "Belowground"), as.factor))
DryMass <- DryMass %>%
dplyr::mutate(across(c("ID", "Plot", "Species", "Aboveground", "Belowground"), as.factor))
sessionInfo()
## R version 4.4.1 (2024-06-14)
## Platform: aarch64-apple-darwin20
## Running under: macOS 15.0.1
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: America/Bogota
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] lme4_1.1-35.5 Matrix_1.7-0 glmmTMB_1.1.9 ggfortify_0.4.17
## [5] survival_3.6-4 ggpubr_0.6.0 ggbeeswarm_0.7.2 report_0.5.8
## [9] reshape2_1.4.4 emmeans_1.10.2 betareg_3.2-0 lmtest_0.9-40
## [13] zoo_1.8-12 pscl_1.5.9 Rmisc_1.5.1 plyr_1.8.9
## [17] lattice_0.22-6 effects_4.2-2 car_3.1-2 carData_3.0-5
## [21] MASS_7.3-60.2 kableExtra_1.4.0 dplyr_1.1.4 ggplot2_3.5.1
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.1 viridisLite_0.4.2 vipor_0.4.7
## [4] fastmap_1.2.0 TH.data_1.1-2 digest_0.6.35
## [7] estimability_1.5.1 lifecycle_1.0.4 magrittr_2.0.3
## [10] compiler_4.4.1 rlang_1.1.4 sass_0.4.9
## [13] tools_4.4.1 utf8_1.2.4 yaml_2.3.8
## [16] knitr_1.47 ggsignif_0.6.4 xml2_1.3.6
## [19] abind_1.4-5 multcomp_1.4-25 numDeriv_2016.8-1.1
## [22] purrr_1.0.2 withr_3.0.0 nnet_7.3-19
## [25] grid_4.4.1 stats4_4.4.1 fansi_1.0.6
## [28] xtable_1.8-4 colorspace_2.1-0 scales_1.3.0
## [31] insight_0.20.1 cli_3.6.3 mvtnorm_1.2-5
## [34] survey_4.4-2 rmarkdown_2.27 generics_0.1.3
## [37] rstudioapi_0.16.0 minqa_1.2.7 DBI_1.2.3
## [40] cachem_1.1.0 stringr_1.5.1 modeltools_0.2-23
## [43] splines_4.4.1 mitools_2.4 vctrs_0.6.5
## [46] boot_1.3-30 sandwich_3.1-0 jsonlite_1.8.8
## [49] rstatix_0.7.2 Formula_1.2-5 beeswarm_0.4.0
## [52] systemfonts_1.1.0 tidyr_1.3.1 jquerylib_0.1.4
## [55] glue_1.8.0 nloptr_2.1.0 codetools_0.2-20
## [58] stringi_1.8.4 gtable_0.3.5 munsell_0.5.1
## [61] tibble_3.2.1 pillar_1.9.0 htmltools_0.5.8.1
## [64] TMB_1.9.14 R6_2.5.1 evaluate_0.24.0
## [67] backports_1.5.0 broom_1.0.6 bslib_0.7.0
## [70] Rcpp_1.0.12 flexmix_2.3-19 gridExtra_2.3
## [73] svglite_2.1.3 coda_0.19-4.1 nlme_3.1-164
## [76] mgcv_1.9-1 xfun_0.45 pkgconfig_2.0.3
search()
## [1] ".GlobalEnv" "package:lme4" "package:Matrix"
## [4] "package:glmmTMB" "package:ggfortify" "package:survival"
## [7] "package:ggpubr" "package:ggbeeswarm" "package:report"
## [10] "package:reshape2" "package:emmeans" "package:betareg"
## [13] "package:lmtest" "package:zoo" "package:pscl"
## [16] "package:Rmisc" "package:plyr" "package:lattice"
## [19] "package:effects" "package:car" "package:carData"
## [22] "package:MASS" "package:kableExtra" "package:dplyr"
## [25] "package:ggplot2" "package:stats" "package:graphics"
## [28] "package:grDevices" "package:utils" "package:datasets"
## [31] "package:methods" "Autoloads" "package:base"
For the WetMass dataframe, we need to still calculate average leaf area. For most plants, only one leaf was measured due to time constraints. We may not be able to use leaf area. First we need to convert the “0” values into “NA”s. Then we can use the rowMeans function to calculate the mean leaf area. Note that the leaf area’s only go up to 16 because leaves IDed above 17 were never measured for area.
WetMass[, c("TotArea.1", "TotArea.2","TotArea.3", "TotArea.4", "TotArea.5","TotArea.6",
"TotArea.7","TotArea.8","TotArea.9","TotArea.10",
"TotArea.11","TotArea.12","TotArea.13","TotArea.14",
"TotArea.15","TotArea.16")] <- lapply(WetMass[, c("TotArea.1",
"TotArea.2","TotArea.3", "TotArea.4", "TotArea.5","TotArea.6",
"TotArea.7","TotArea.8","TotArea.9","TotArea.10",
"TotArea.11","TotArea.12","TotArea.13","TotArea.14",
"TotArea.15","TotArea.16")], function(x) ifelse(x == 0, NA, x))
WetMass$AveArea <- rowMeans(WetMass[, c("TotArea.1", "TotArea.2",
"TotArea.3", "TotArea.4", "TotArea.5","TotArea.6",
"TotArea.7","TotArea.8","TotArea.9","TotArea.10",
"TotArea.11","TotArea.12","TotArea.13","TotArea.14",
"TotArea.15","TotArea.16")], na.rm = TRUE)
We will to fully understand our variables before we launch into
statistical analysis. Let’s get somebasic descriptive statisticsgoing
with the summary function.
summary(Fern)
## ID Plot Species Aboveground
## 1 : 2 Shadehouse:85 Cymbopetalum baillonii:294 Covered :240
## 10 : 2 15.3 :80 Sapindus saponaria :271 Exposed :240
## 100 : 2 15.4 :80 Invernadero: 85
## 101 : 2 24.2 :80
## 102 : 2 6.3 :80
## 103 : 2 8.2 :80
## (Other):553 (Other) :80
## Belowground Germinated Germination.Date Death.Date
## Invernadero: 85 No :406 Min. :2019-07-12 Min. :2019-07-19
## Protected :240 Yes:159 1st Qu.:2019-07-19 1st Qu.:2019-08-31
## Unprotected:240 Median :2019-07-26 Median :2019-11-29
## Mean :2019-09-12 Mean :2020-04-17
## 3rd Qu.:2019-08-16 3rd Qu.:2021-04-16
## Max. :2021-03-26 Max. :2021-05-25
## NA's :407 NA's :463
## Status Weeks.Alive Coty.Date Weeks.Until.Coty
## Min. :0.0000 Min. : 1.00 Min. :2019-07-19 Min. :1.00
## 1st Qu.:0.0000 1st Qu.: 8.75 1st Qu.:2019-08-02 1st Qu.:1.00
## Median :1.0000 Median :51.00 Median :2019-08-16 Median :2.00
## Mean :0.6497 Mean :50.81 Mean :2019-10-20 Mean :2.35
## 3rd Qu.:1.0000 3rd Qu.:93.00 3rd Qu.:2019-09-20 3rd Qu.:3.00
## Max. :1.0000 Max. :98.00 Max. :2021-04-02 Max. :7.00
## NA's :408 NA's :409 NA's :462 NA's :462
## Date.First.Leaf Date.Highest.Leaf Weeks.Between.Leaves
## Min. :2019-07-19 Min. :2019-08-16 Min. : 0.00
## 1st Qu.:2019-08-16 1st Qu.:2019-11-23 1st Qu.: 5.25
## Median :2019-08-23 Median :2020-10-02 Median :36.00
## Mean :2019-10-31 Mean :2020-08-06 Mean :40.02
## 3rd Qu.:2019-10-26 3rd Qu.:2021-04-09 3rd Qu.:76.50
## Max. :2021-04-09 Max. :2021-05-21 Max. :94.00
## NA's :435 NA's :435 NA's :435
## LeafRate State StartLeaves MaxLeaves
## Min. :0.0500 Dead : 80 Min. :0.0000 Min. : 0.000
## 1st Qu.:0.1925 Dead Harvest: 24 1st Qu.:0.0000 1st Qu.: 0.000
## Median :0.2950 Harvested : 55 Median :0.0000 Median : 0.000
## Mean :0.3887 MIA :406 Mean :0.4212 Mean : 2.701
## 3rd Qu.:0.4675 3rd Qu.:0.0000 3rd Qu.: 0.000
## Max. :2.3300 Max. :4.0000 Max. :44.000
## NA's :459
## LastLeaves
## Min. : 2.00
## 1st Qu.: 5.00
## Median : 8.00
## Mean :11.62
## 3rd Qu.:16.00
## Max. :40.00
## NA's :510
Variable description:
ID = ID of plant per species.Plot = Actually subplot, location where seeds were
sowed.Species = Species planted (two levels)Aboveground = Experimental treatment (two levels) +
ShadehouseBelowground = Experimental treatment (two levels) +
ShadehouseGerminated = Boolean.Germination.Date = Date when germinating seed was first
observed.Death.Date = Date when plant was not considered alive
anymore.Weeks.Alive = Numbers of weeks between
Germination.Date and Death.Date.Coty.Date = Date when cotyledons were first
observed.Weeks.Until.Coty = Numbers of weeks between
Germination.Date and Coty.Date.Date.First.Leaf = Date when the first leaf was
observed.Date.Highest.Leaf = Date when the highest leaf count
was observed.Weeks.Between.Leaves = Numbers of weeks between
Date.First.Leaf and Date.Highest.Leaf.LeafRate = Rate of leaf production calculated using
Weeks.Between.Leaves and
Date.Highest.Leaf.State = Four level factor. MIA indicates
the plant never germinated. Dead Harvest indicates that the
plant was harvested but was dead by then. Harvested and
Dead self-explanatory.StartLeaves = Number of leaves when leaves first
reported.MaxLeaves = Highest number of leaves recorded on
plant.LastLeaves = Number of leaves when last alive.Let’s make a TRUE ID variable.
Fern$TrueID <- paste(Fern$ID, Fern$Species)
Fern$Treatment <- as.factor(paste(Fern$Aboveground, "+", Fern$Belowground))
WetMass$TrueID <- paste(WetMass$ID, WetMass$Species)
DryMass$TrueID <- paste(DryMass$ID, DryMass$Species)
How many plots do we have? Which ones?
T1 <- as.data.frame(as.table(summary(Fern$Plot)))
colnames(T1) <- c("Subplot", "Count")
kbl(T1, format = "html",
table.attr = "style='width:35%;'") %>%
kable_classic()
| Subplot | Count |
|---|---|
| 10.2 | 60 |
| 10.3 | 20 |
| 15.3 | 80 |
| 15.4 | 80 |
| 24.2 | 80 |
| 6.3 | 80 |
| 8.2 | 80 |
| Shadehouse | 85 |
Answer: We have six plots in our study: Plot 6, Plot 8, Plot 10, Plot 15, and Plot 24. Regarding Plot 10, two different subplots were used to achieve the full sample size of 80. In plot 15 we had two replicates. We also have a Shadehouse category.
levels(Fern$Plot) <- list("10" = "10.2", "10" = "10.3", "15A" = "15.3",
"15B" = "15.4", "24" = "24.2", "6" = "6.3",
"8" = "8.2", "Shadehouse" = "Shadehouse")
levels(Fern$Plot)
## [1] "10" "15A" "15B" "24" "6"
## [6] "8" "Shadehouse"
Let’s create a dataframe that excludes the Shadehouse and another one
that only holds Shadehouse. We’ll remove them using the combo of
subset and droplevels. We’ll also create
dataframes for each species.
FernExp <- droplevels(subset(Fern, Fern$Plot != "Shadehouse"))
Shadehouse <- droplevels(subset(Fern, Fern$Plot == "Shadehouse"))
FernExp.SASA <- droplevels(subset(FernExp, FernExp$Species == "Sapindus saponaria"))
FernExp.CYBA <- droplevels(subset(FernExp, FernExp$Species == "Cymbopetalum baillonii"))
FernExp.SASA$Belowground <- relevel(FernExp.SASA$Belowground, ref = "Unprotected")
FernExp.CYBA$Belowground <- relevel(FernExp.CYBA$Belowground, ref = "Unprotected")
How many species did we plant? Which ones and in what treatments?
T2 <- as.data.frame(table(FernExp$Species, FernExp$Aboveground, FernExp$Belowground))
colnames(T2) <- c("Species", "Aboveground", "Belowground", "Count")
kbl(T2, format = "html",
table.attr = "style='width:75%;'") %>%
kable_classic()
| Species | Aboveground | Belowground | Count |
|---|---|---|---|
| Cymbopetalum baillonii | Covered | Protected | 60 |
| Sapindus saponaria | Covered | Protected | 60 |
| Cymbopetalum baillonii | Exposed | Protected | 60 |
| Sapindus saponaria | Exposed | Protected | 60 |
| Cymbopetalum baillonii | Covered | Unprotected | 60 |
| Sapindus saponaria | Covered | Unprotected | 60 |
| Cymbopetalum baillonii | Exposed | Unprotected | 60 |
| Sapindus saponaria | Exposed | Unprotected | 60 |
Interpretation: We planted two species, Cymbopetalum baillonii and Sapindus saponaria.
The experiment follows a 2x2 Factorial Design using two variables: Aboveground (Covered vs Exposed) and Belowground (Protected vs Unprotected). Seedlings were planted in plots that had low tree canopy cover and high fern density (see Beltran et al. 2020). Natural conditions for seedlings growing in these conditions would involve a combination of the Covered and Unprotected levels.
In this experiment, Exposed was accomplished by tying up the covering fronds in nylon and pulling them away. Protected was accomplished by germinating seedlings in 15 cm diameter PVC tubes with a 1 mm mosquito netting attached to the bottom via rubber bands (height of tube 20 cm).
Before we get started, it is important to note that in some cases, plants burst through the netting that composed the treatment that limited belowground competition. In some cases, the fern the tubes from above or below. Let’s take a look at some descriptive statistics to decide whether these should be excluded.
exclude_ids <- c(83, 84, 87, 88, 167)
WetMass2 <- WetMass %>%
filter(!(ID %in% exclude_ids))
DryMass2 <- DryMass %>%
filter(!(ID %in% exclude_ids))
WetMassExp <- droplevels(subset(WetMass2, WetMass$Plot != "Invernadero"))
DryMassExp <- droplevels(subset(DryMass2, DryMass$Plot != "Invernadero"))
WetMassExp$Treatment <- as.factor(paste(WetMassExp$Aboveground, "+", WetMassExp$Belowground))
DryMassExp$Treatment <- as.factor(paste(DryMassExp$Aboveground, "+", DryMassExp$Belowground))
levels(WetMassExp$Treatment) <- list("No Competition" = "Exposed + Protected",
"Belowground Competition" = "Exposed + Unprotected",
"Aboveground Competition" = "Covered + Protected",
"Full Competition" = "Covered + Unprotected")
levels(DryMassExp$Treatment) <- list("No Competition" = "Exposed + Protected",
"Belowground Competition" = "Exposed + Unprotected",
"Aboveground Competition" = "Covered + Protected",
"Full Competition" = "Covered + Unprotected")
WetMassExp$Harvest.Notes <- as.factor(WetMassExp$Harvest.Notes)
summarySE(WetMassExp, measurevar="TotalWetMass", groupvars=c("Species", "Harvest.Notes", "Treatment"),
na.rm = TRUE)[c(3, 5, 6),-c(6,8)]
## Species Harvest.Notes Treatment N
## 3 Cymbopetalum baillonii Aboveground Competition 4
## 5 Cymbopetalum baillonii Died / Drying <NA> 1
## 6 Cymbopetalum baillonii Fern burst through Aboveground Competition 2
## TotalWetMass se
## 3 1.9225 0.7297645
## 5 12.3500 NA
## 6 1.9850 0.3650000
summarySE(WetMassExp, measurevar="TotalWetMass", groupvars=c("Species", "Harvest.Notes", "Treatment"),
na.rm = TRUE)[c(1,7),-c(6,8)]
## Species Harvest.Notes Treatment N TotalWetMass
## 1 Cymbopetalum baillonii No Competition 3 6.72
## 7 Cymbopetalum baillonii Fern invaded Aboveground Competition 1 2.22
## se
## 1 3.821941
## 7 NA
summarySE(WetMassExp, measurevar="TotalWetMass", groupvars=c("Species", "Harvest.Notes", "Treatment"),
na.rm = TRUE)[c(10, 12),-c(6,8)]
## Species Harvest.Notes Treatment N TotalWetMass
## 10 Sapindus saponaria Aboveground Competition 7 0.8357143
## 12 Sapindus saponaria Died / Drying <NA> 4 10.2500000
## se
## 10 0.1581418
## 12 3.1677568
summarySE(WetMassExp, measurevar="TotalWetMass", groupvars=c("Species", "Harvest.Notes", "Treatment"),
na.rm = TRUE)[c(8, 13),-c(6,8)]
## Species Harvest.Notes Treatment N TotalWetMass
## 8 Sapindus saponaria No Competition 9 22.17111
## 13 Sapindus saponaria Fern burst through Aboveground Competition 1 3.32000
## se
## 8 9.670179
## 13 NA
Interpretation:
Fresh biomass of Cymbopetalum in the aboveground only treatment was not affected by the invading ferns into the tubes.
Fresh biomass of Cymbopetalum in the no competition treatment differed greatly when the plants broke through the tubes.
Fresh biomass of Sapindus saponaria in the aboveground only treatment differed, but only one plant was affected by the fern invasion, and oddly had more mass.
Fresh biomass of Sapindus saponaria in the competition treatment was larger when the plant broke through.
Conclusion: Perhaps we should exclude from the biomass and mortality analyses the plants that broke the netting treatment. IDs: 83, 84, 87, 88, and 167.
How many seeds germinated in the experimental plots?
T3 <- as.data.frame(table(FernExp$Germinated))
colnames(T3) <- c("Germinated", "Frequency")
T3$Percent <- round(T3$Frequency / 480 * 100, 2)
kbl(T3, format = "html",
table.attr = "style='width:40%;'") %>% kable_classic()
| Germinated | Frequency | Percent |
|---|---|---|
| No | 350 | 72.92 |
| Yes | 130 | 27.08 |
Answer: 130 out of 480 seeds sowed in the experimental plots germinated (27.08%)
How many seeds germinated per species in the experimental plots?
T4 <- as.data.frame(table(FernExp$Germinated, FernExp$Species))
colnames(T4) <- c("Germinated", "Species", "Frequency")
T4$Percent <- round(T4$Frequency / 240 * 100, 2)
T4 <- droplevels(subset(T4, T4$Germinated == "Yes"))
kbl(T4[,-1], format = "html",
table.attr = "style='width:50%;'", row.names = FALSE) %>% kable_classic()
| Species | Frequency | Percent |
|---|---|---|
| Cymbopetalum baillonii | 82 | 34.17 |
| Sapindus saponaria | 48 | 20.00 |
Answer: A total of 82 Cymbopetalum bailloni seeds germinated (34.17%) vs. 48 Sapindus saponaria seeds (20%).
How many seeds germinated per treatment and species?
Lets create a new variable that reflects the combination of treatments.
FernExp$Treatment <- as.factor(paste(FernExp$Aboveground, "+", FernExp$Belowground))
levels(FernExp$Treatment) <- list("No Competition" = "Exposed + Protected",
"Belowground Competition" = "Exposed + Unprotected",
"Aboveground Competition" = "Covered + Protected",
"Full Competition" = "Covered + Unprotected")
T5 <- as.data.frame(table(FernExp$Germinated, FernExp$Species, FernExp$Treatment))
colnames(T5) <- c("Germinated", "Species", "Treatment", "Frequency")
T5$Percent <- round(T5$Frequency / 60 ,4) *100
T5 <- droplevels(subset(T5, T5$Germinated == "Yes"))
kbl(T5[,-1], format = "html",
table.attr = "style='width:55%;'", row.names = FALSE) %>% kable_classic()
| Species | Treatment | Frequency | Percent |
|---|---|---|---|
| Cymbopetalum baillonii | No Competition | 17 | 28.33 |
| Sapindus saponaria | No Competition | 15 | 25.00 |
| Cymbopetalum baillonii | Belowground Competition | 12 | 20.00 |
| Sapindus saponaria | Belowground Competition | 11 | 18.33 |
| Cymbopetalum baillonii | Aboveground Competition | 22 | 36.67 |
| Sapindus saponaria | Aboveground Competition | 10 | 16.67 |
| Cymbopetalum baillonii | Full Competition | 31 | 51.67 |
| Sapindus saponaria | Full Competition | 12 | 20.00 |
How many seeds germinated species in the Shadehouse?
T6 <- as.data.frame(table(Shadehouse$Germinated, Shadehouse$Species))
colnames(T6) <- c("Germinated", "Species", "Count")
T6$Total <- as.integer(c("54", "54", "31", "31"))
T6$Percent <- round(T6$Count / T6$Total,4) * 100
T6 <- droplevels(subset(T6, T6$Germinated == "Yes"))
kbl(T6[,-c(1,4)], format = "html",
table.attr = "style='width:50%;'", row.names = FALSE) %>% kable_classic()
| Species | Count | Percent |
|---|---|---|
| Cymbopetalum baillonii | 25 | 46.3 |
| Sapindus saponaria | 4 | 12.9 |
Answer: Highest and lowest germination in the
experimental treatments for Cymbopetalum was in the
Covered + Unprotected (25.83%) and
Exposed + Unprotected (10%) treatments respectively.
Shadehouse germination was highest at 46.3%. For S. saponaria,
highest germination was achieved under the
Exposed + Protected (12.5%) and lowest under the
Covered + Protected (8.33%) treatment. Shadehouse
germination was highest at 12.9%.
Let’s get some averages and standard error.
T7 <- as.data.frame(table(FernExp$Germinated, FernExp$Plot, FernExp$Treatment, FernExp$Species))
T7$Percent <- T7$Freq/10 *100
colnames(T7) <- c("Germinated", "Plot", "Treatment", "Species", "Count", "Percent")
T7 <- droplevels(subset(T7, T7$Germinated == "Yes"))
kbl(summarySE(T7, measurevar="Percent", groupvars=c("Species", "Treatment"),
na.rm = TRUE)[,-c(3,5,7)], format = "html",
table.attr = "style='width:75%;'") %>% kable_classic()
| Species | Treatment | Percent | se |
|---|---|---|---|
| Cymbopetalum baillonii | No Competition | 28.33333 | 9.098229 |
| Cymbopetalum baillonii | Belowground Competition | 20.00000 | 9.309493 |
| Cymbopetalum baillonii | Aboveground Competition | 36.66667 | 6.666667 |
| Cymbopetalum baillonii | Full Competition | 51.66667 | 8.333333 |
| Sapindus saponaria | No Competition | 25.00000 | 4.281744 |
| Sapindus saponaria | Belowground Competition | 18.33333 | 4.772607 |
| Sapindus saponaria | Aboveground Competition | 16.66667 | 4.944132 |
| Sapindus saponaria | Full Competition | 20.00000 | 6.831301 |
T7B <- as.data.frame(table(FernExp$Germinated, FernExp$Plot,
FernExp$Aboveground, FernExp$Belowground, FernExp$Species))
colnames(T7B) <- c("Germinated", "Plot", "Aboveground", "Belowground", "Species", "Count")
T7B <- droplevels(subset(T7B, T7B$Germinated == "Yes"))
T7B$TRUEID <- paste(T7B$Plot, T7B$Aboveground, T7B$Belowground, T7B$Species)
Answer: Cymbopetalum showed highest mean
germination rates under Full Competition, followed by
Aboveground Competition, No Competition, and
lastly Belowground Competition. Sapindus showed
highest mean germination rates under No Competition,
followed by Full Competition, then
Belowground Competition, and lastly
Aboveground Competition.
How many seeds that germinated in the experimental plots died?
T8 <- as.data.frame(table(FernExp$State))
colnames(T8) <- c("State", "Died")
T8$Percent <- round(T8$Died / 130 * 100, 2)
T8 <- droplevels(subset(T8, T8$State == "Dead"))
kbl(T8, format = "html",
table.attr = "style='width:30%;'") %>% kable_classic()
| State | Died | Percent |
|---|---|---|
| Dead | 75 | 57.69 |
Answer: 75 out of 130 seedlings that germinated in the experimental plots died within the growth period (57.69%)
How many seedlings per species died in the experimental plots?
T9 <- as.data.frame(table(FernExp$State, FernExp$Species))
colnames(T9) <- c("State", "Species", "Died")
T9 <- droplevels(subset(T9, T9$State == "Dead"))
T9$Percent <- round(T9$Died / T4$Frequency, 4) * 100
kbl(T9[,-1], format = "html",
table.attr = "style='width:50%;'", row.names = FALSE) %>% kable_classic()
| Species | Died | Percent |
|---|---|---|
| Cymbopetalum baillonii | 59 | 71.95 |
| Sapindus saponaria | 16 | 33.33 |
32/55*100
## [1] 58.18182
Answer: A total of 59 Cymbopetalum bailloni seedlings died out of 82 (71.95%) vs. 16 Sapindus saponaria seedlings out of 48 (33.33%).
How many seeds died per treatment? How many in the Shadehouse?
T10 <- as.data.frame(table(FernExp$State, FernExp$Species, FernExp$Treatment))
colnames(T10) <- c("State", "Species", "Treatment", "Died")
T10 <- droplevels(subset(T10, T10$State == "Dead"))
T10$Percent <- round(T10$Died / T5$Frequency ,4) *100
kbl(T10[,-1], format = "html",
table.attr = "style='width:55%;'", row.names = FALSE) %>% kable_classic()
| Species | Treatment | Died | Percent |
|---|---|---|---|
| Cymbopetalum baillonii | No Competition | 10 | 58.82 |
| Sapindus saponaria | No Competition | 6 | 40.00 |
| Cymbopetalum baillonii | Belowground Competition | 9 | 75.00 |
| Sapindus saponaria | Belowground Competition | 4 | 36.36 |
| Cymbopetalum baillonii | Aboveground Competition | 14 | 63.64 |
| Sapindus saponaria | Aboveground Competition | 2 | 20.00 |
| Cymbopetalum baillonii | Full Competition | 26 | 83.87 |
| Sapindus saponaria | Full Competition | 4 | 33.33 |
T11 <- as.data.frame(table(Shadehouse$State, Shadehouse$Species))
colnames(T11) <- c("State", "Species", "Died")
T11 <- droplevels(subset(T11, T11$State != "MIA"))
T11$Total <- as.numeric(c("25", "25", "4", "4"))
T11$Percent <- round(T11$Died / T11$Total,4) * 100
kbl(T11[,-c(4)], format = "html",
table.attr = "style='width:50%;'", row.names = FALSE) %>% kable_classic()
| State | Species | Died | Percent |
|---|---|---|---|
| Dead | Cymbopetalum baillonii | 5 | 20 |
| Dead Harvest | Cymbopetalum baillonii | 20 | 80 |
| Dead | Sapindus saponaria | 0 | 0 |
| Dead Harvest | Sapindus saponaria | 4 | 100 |
Answer: Hmm.. well, its complicated. We had a problem on the last month of the experiment. Unbeknownst to me, the field station workers stopped watering my plants. My assistant knew, but thought it not important to tell me. By the time I rolled around for the harvest, they were all dead. Those listed as “Dead Harvest” are those that appear to have died as a direct cause of the lack of water and we harvested them. However, this data is hardly comparable, as the leaf data indicates that they had been dying for weeks.
What do mean values of mortality look like?
T12 <- as.data.frame(table(FernExp$State, FernExp$Plot, FernExp$Treatment, FernExp$Species))
colnames(T12) <- c("State", "Plot", "Treatment", "Species", "Count")
T12 <- droplevels(subset(T12, T12$State == "Dead"))
T12$TRUEID <- paste(T12$Plot, T12$Treatment, T12$Species)
T7$TRUEID <- paste(T7$Plot, T7$Treatment, T7$Species)
T12$Germinated <- with(T7, Count[match(T12$TRUEID, TRUEID)])
T12 <- droplevels(subset(T12, T12$Germinated != "0"))
T12$ProportionDead <- round(T12$Count / T12$Germinated,2)
kbl(summarySE(T12, measurevar="ProportionDead", groupvars=c("Species", "Treatment"),
na.rm = TRUE)[,-c(5,7)], format = "html",
table.attr = "style='width:75%;'") %>% kable_classic()
| Species | Treatment | N | ProportionDead | se |
|---|---|---|---|---|
| Cymbopetalum baillonii | No Competition | 5 | 0.6340000 | 0.1611087 |
| Cymbopetalum baillonii | Belowground Competition | 4 | 0.7925000 | 0.1247247 |
| Cymbopetalum baillonii | Aboveground Competition | 6 | 0.6366667 | 0.1569855 |
| Cymbopetalum baillonii | Full Competition | 6 | 0.8883333 | 0.0715736 |
| Sapindus saponaria | No Competition | 6 | 0.4166667 | 0.1598054 |
| Sapindus saponaria | Belowground Competition | 6 | 0.2500000 | 0.1118034 |
| Sapindus saponaria | Aboveground Competition | 5 | 0.1320000 | 0.0808332 |
| Sapindus saponaria | Full Competition | 6 | 0.4216667 | 0.1899196 |
T12B <- as.data.frame(table(FernExp$State, FernExp$Plot,
FernExp$Aboveground, FernExp$Belowground, FernExp$Species))
colnames(T12B) <- c("State", "Plot", "Aboveground", "Belowground", "Species", "Count")
T12B <- droplevels(subset(T12B, T12B$State == "Dead"))
T12B$TRUEID <- paste(T12B$Plot, T12B$Aboveground, T12B$Belowground, T12B$Species)
T12B$Germinated <- with(T7B, Count[match(T12B$TRUEID, TRUEID)])
T12B <- droplevels(subset(T12B, T12B$Germinated != "0"))
T12B$ProportionDead <- round(T12B$Count / T12B$Germinated,2)
Answer: Mortality for Cymbopetalum was
highest under Full Competition, followed by
Aboveground Competition, far from
No Competition and Belowground Competition.
Mortality for Sapindus was highest under
Belowground Competition, followed by
Aboveground Competition and No Competition,
lastly Full Competition.
Question: How long did each species live on average?
kbl(summarySE(FernExp, measurevar="Weeks.Alive", groupvars=c("Species"), na.rm = TRUE)[,-c(4,6)]) %>% kable_classic()
| Species | N | Weeks.Alive | se |
|---|---|---|---|
| Cymbopetalum baillonii | 81 | 37.92593 | 4.346411 |
| Sapindus saponaria | 46 | 55.30435 | 4.990994 |
Answer: On average, C. baillonii seedlings lived 37.93 weeks out of 97 possible weeks, while S. saponaria lived 55.3 weeks.
On average, how long did seedlings lived by species and treatment?
kbl(summarySE(FernExp, measurevar="Weeks.Alive", groupvars=c("Species", "Aboveground","Belowground"),
na.rm = TRUE)[,-c(4,6,8)], format = "html",
table.attr = "style='width:75%;'") %>% kable_classic()
| Species | Aboveground | Belowground | Weeks.Alive | se |
|---|---|---|---|---|
| Cymbopetalum baillonii | Covered | Protected | 46.72727 | 9.085496 |
| Cymbopetalum baillonii | Covered | Unprotected | 31.50000 | 6.367978 |
| Cymbopetalum baillonii | Exposed | Protected | 43.41176 | 10.826526 |
| Cymbopetalum baillonii | Exposed | Unprotected | 30.08333 | 9.895253 |
| Sapindus saponaria | Covered | Protected | 58.20000 | 8.807572 |
| Sapindus saponaria | Covered | Unprotected | 41.91667 | 9.227313 |
| Sapindus saponaria | Exposed | Protected | 59.00000 | 10.327601 |
| Sapindus saponaria | Exposed | Unprotected | 63.30000 | 10.987923 |
Answer: Seedlings that lived the longest belonged to
the Protected treatment. Full Competition showed the lowest
values.
What does the distribution of the number of weeks that seedlings lived look like?
ggplot(FernExp, aes(x=Weeks.Alive)) +
geom_histogram(aes(y=..density..), colour="black", fill="cornsilk") +
geom_density(alpha=.2, fill="#FF6666") + # easier to see distribution
geom_vline(aes(xintercept=mean(Weeks.Alive, na.rm=TRUE)),color="cadetblue4",
linetype="dashed", size=1) + # This adds a blue line where the mean is.
xlab("Weeks Alive") +
ylab("Density") +
theme_bw(base_size = 15)
Answer: We do not have a clear distribution either… we’ll have to take a look at our classic model shapes.
Did treatment affect the max number of leaves? We’ll need to exclude only those that germinated.
FernExp.G <- droplevels(subset(FernExp, FernExp$Germinated == "Yes"))
Let’s exclude the plants that burst through the netting. IDs: 83, 84, 87, 88, and 167.
FernExp.G2 <- FernExp.G %>%
filter(!(ID %in% exclude_ids))
FernExp.G.CYBA <- droplevels(subset(FernExp.G2, FernExp.G2$Species == "Cymbopetalum baillonii"))
FernExp.G.SASA <- droplevels(subset(FernExp.G2, FernExp.G2$Species == "Sapindus saponaria"))
options(digits = 5)
kbl(summarySE(FernExp.G2, measurevar="MaxLeaves", groupvars=c("Species", "Treatment"), na.rm = TRUE)[,-c(3,5,7)]) %>% kable_classic()
| Species | Treatment | MaxLeaves | se |
|---|---|---|---|
| Cymbopetalum baillonii | No Competition | 4.7692 | 2.34563 |
| Cymbopetalum baillonii | Belowground Competition | 4.6667 | 1.34465 |
| Cymbopetalum baillonii | Aboveground Competition | 5.3182 | 1.12313 |
| Cymbopetalum baillonii | Full Competition | 3.5484 | 0.87109 |
| Sapindus saponaria | No Competition | 10.0714 | 2.64464 |
| Sapindus saponaria | Belowground Competition | 10.8182 | 2.76609 |
| Sapindus saponaria | Aboveground Competition | 5.3000 | 0.55877 |
| Sapindus saponaria | Full Competition | 5.5833 | 1.31690 |
Answer: For C. baillonii, the lowest max
leaf count was 3.55 ± 0.871 when Full Competition and
highest with 10.07 ± 2.644 when No Competition. For S.
saponaria, lowest max leaf count was 5.30 ± 0.559 when
Full Competition and highest max leaf count was 10.82 ±
2.766 when Belowground Competition, followed closely by
No Competition 10.07 ± 2.645.
What about in the shadehouse?
Shadehouse.G <- droplevels(subset(Shadehouse, Shadehouse$Germinated == "Yes"))
summarySE(Shadehouse.G, measurevar="MaxLeaves", groupvars=c("Species"), na.rm = TRUE)[,-c(2,4,6)]
## Species MaxLeaves se
## 1 Cymbopetalum baillonii 23.2 1.8956
## 2 Sapindus saponaria 11.5 1.6583
Did seedling height differ between species and treatment?
kbl(summarySE(WetMassExp, measurevar="Height..cm.", groupvars=c("Aboveground"),
na.rm = TRUE)[,-c(4,6)], format = "html",
table.attr = "style='width:75%;'") %>% kable_classic()
| Aboveground | N | Height..cm. | se |
|---|---|---|---|
| Covered | 28 | 16.536 | 1.1768 |
| Exposed | 22 | 29.045 | 5.7545 |
| Invernadero | 0 | NaN | NA |
kbl(summarySE(WetMassExp, measurevar="Height..cm.", groupvars=c("Aboveground", "Species"),
na.rm = TRUE)[,-c(3,6)], format = "html",
table.attr = "style='width:75%;'") %>% kable_classic()
| Aboveground | Species | Height..cm. | sd | ci |
|---|---|---|---|---|
| Covered | Cymbopetalum baillonii | 19.00 | 5.7300 | 3.4626 |
| Covered | Sapindus saponaria | 14.40 | 6.0095 | 3.3280 |
| Exposed | Cymbopetalum baillonii | 20.50 | 10.5024 | 11.0216 |
| Exposed | Sapindus saponaria | 32.25 | 30.7083 | 16.3633 |
| Invernadero | Cymbopetalum baillonii | NaN | NA | NA |
| Invernadero | Sapindus saponaria | NaN | NA | NA |
14.40/32.25
## [1] 0.44651
WetMassExp <- WetMassExp %>%
filter(!is.na(Treatment))
kbl(summarySE(WetMassExp, measurevar="Height..cm.", groupvars=c("Species", "Treatment"),
na.rm = TRUE)[,-c(5,7)], format = "html",
table.attr = "style='width:75%;'") %>% kable_classic()
| Species | Treatment | N | Height..cm. | se |
|---|---|---|---|---|
| Cymbopetalum baillonii | No Competition | 3 | 19.333 | 5.8119 |
| Cymbopetalum baillonii | Belowground Competition | 3 | 21.667 | 7.5351 |
| Cymbopetalum baillonii | Aboveground Competition | 7 | 18.429 | 2.1476 |
| Cymbopetalum baillonii | Full Competition | 6 | 19.667 | 2.5517 |
| Sapindus saponaria | No Competition | 9 | 29.889 | 6.7401 |
| Sapindus saponaria | Belowground Competition | 7 | 35.286 | 16.0055 |
| Sapindus saponaria | Aboveground Competition | 8 | 13.000 | 2.1381 |
| Sapindus saponaria | Full Competition | 7 | 16.000 | 2.2678 |
Should we be excluding seedlings that germinated late? Not a huge difference in numbers here.
summary(FernExp$Germination.Date)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## "2019-07-12" "2019-07-19" "2019-08-02" "2019-09-21" "2019-09-06" "2021-03-26"
## NA's
## "351"
WetMassExp$Weeks.Alive <- with(FernExp, Weeks.Alive[match(WetMassExp$TrueID, TrueID)])
WetMassExp$Germination.Date <- with(FernExp, Germination.Date[match(WetMassExp$TrueID, TrueID)])
WetMassExp.B <- WetMassExp
WetMassExp.B <- droplevels(subset(WetMassExp.B, WetMassExp.B$Germination.Date != "2021-03-26"))
WetMassExp.B <- droplevels(subset(WetMassExp.B, WetMassExp.B$Germination.Date != "2020-07-17"))
WetMassExp.B <- droplevels(subset(WetMassExp.B, WetMassExp.B$Germination.Date != "2020-07-03"))
WetMassExp.B <- droplevels(subset(WetMassExp.B, WetMassExp.B$Germination.Date != "2020-06-05"))
kbl(summarySE(WetMassExp.B, measurevar="Height..cm.", groupvars=c("Species", "Treatment"),
na.rm = TRUE)[,-c(5,7)], format = "html",
table.attr = "style='width:75%;'") %>% kable_classic()
| Species | Treatment | N | Height..cm. | se |
|---|---|---|---|---|
| Cymbopetalum baillonii | No Competition | 3 | 19.333 | 5.8119 |
| Cymbopetalum baillonii | Belowground Competition | 3 | 21.667 | 7.5351 |
| Cymbopetalum baillonii | Aboveground Competition | 7 | 18.429 | 2.1476 |
| Cymbopetalum baillonii | Full Competition | 5 | 20.200 | 3.0561 |
| Sapindus saponaria | No Competition | 8 | 27.375 | 7.0911 |
| Sapindus saponaria | Belowground Competition | 7 | 35.286 | 16.0055 |
| Sapindus saponaria | Aboveground Competition | 7 | 13.429 | 2.4188 |
| Sapindus saponaria | Full Competition | 5 | 18.000 | 2.3875 |
Answer: Not a huge difference in numbers here. Lets keep them all to maintain our sample size.
Stem Diameter differ between
Aboveground Treatments, Belowground Treatments
and Species?kbl(summarySE(WetMassExp, measurevar="Diameter.Base..cm.", groupvars=c("Species", "Treatment"),
na.rm = TRUE)[,-c(5,7)], format = "html",
table.attr = "style='width:75%;'") %>% kable_classic()
| Species | Treatment | N | Diameter.Base..cm. | se |
|---|---|---|---|---|
| Cymbopetalum baillonii | No Competition | 3 | 0.43333 | 0.08819 |
| Cymbopetalum baillonii | Belowground Competition | 3 | 0.66667 | 0.27285 |
| Cymbopetalum baillonii | Aboveground Competition | 7 | 0.35714 | 0.02020 |
| Cymbopetalum baillonii | Full Competition | 6 | 0.35000 | 0.04282 |
| Sapindus saponaria | No Competition | 9 | 0.72222 | 0.10512 |
| Sapindus saponaria | Belowground Competition | 7 | 2.24286 | 1.33984 |
| Sapindus saponaria | Aboveground Competition | 8 | 0.28750 | 0.02950 |
| Sapindus saponaria | Full Competition | 7 | 0.32857 | 0.03595 |
Interpretation: Yes. For C. baillonii, the
thicker diameters were found in the
Only Belowground competition group (0.67 ± 0.27), followed
by No Competition (0.433 ± 0.088).
Only Aboveground Competition and
Full Competition were about the same (0.357 ± 0.02 and 0.35
± 0.043, respectively.)
For S. saponaria, thickest diameters were found in the
Only Belowground Competition group (2.242 ± 1.339),
followed by No Competition (0.722 ± 0.105),
Full Competition (0.329 ± 0.036), and
Only Aboveground Competition (0.288 ± 0.029).
From harvested plants in the experiment, what was the average wet plant mass by species and treatment?
kbl(summarySE(WetMassExp, measurevar="TotalWetMass", groupvars=c("Species", "Treatment"),
na.rm = TRUE)[,-c(5,7)], format = "html",
table.attr = "style='width:75%;'") %>% kable_classic()
| Species | Treatment | N | TotalWetMass | se |
|---|---|---|---|---|
| Cymbopetalum baillonii | No Competition | 3 | 6.7200 | 3.82194 |
| Cymbopetalum baillonii | Belowground Competition | 3 | 8.3933 | 5.85457 |
| Cymbopetalum baillonii | Aboveground Competition | 7 | 1.9829 | 0.40024 |
| Cymbopetalum baillonii | Full Competition | 6 | 10.1567 | 7.57114 |
| Sapindus saponaria | No Competition | 9 | 22.1711 | 9.67018 |
| Sapindus saponaria | Belowground Competition | 7 | 58.9971 | 46.53081 |
| Sapindus saponaria | Aboveground Competition | 8 | 1.1462 | 0.33940 |
| Sapindus saponaria | Full Competition | 7 | 1.9214 | 0.49068 |
Answer: For C. baillonii, the
Full Competition brought the greatest plant mass (10.157 ±
7.571), followed by Belowground Competition (8.393 ±
5.855), No Competition (6.720 ± 3.822), and lastly
Only Aboveground Competition (1.983 ± 0.40).
For S. saponaria, Only Belowground Competition
brought the greatest plant mass, but with incredible variation (59.00 ±
46.53), followed by No Competition (22.17 ± 9.67).
Full Competition and Aboveground Competition
yielded very low values (1.921 ± 0.49; 1.146 ± 0.339, respectively).
From harvested plants in the experiment, what was the average dry plant mass by species and treatment?
kbl(summarySE(DryMassExp, measurevar="TotalDryMass", groupvars=c("Aboveground"),
na.rm = TRUE)[,-c(4,6)], format = "html",
table.attr = "style='width:75%;'") %>% kable_classic()
| Aboveground | N | TotalDryMass | se |
|---|---|---|---|
| Covered | 28 | 0.79679 | 0.11523 |
| Exposed | 22 | 17.39273 | 9.49590 |
| Invernadero | 5 | 7.13400 | 2.21215 |
(17.392727 - 0.796786) / 0.796786 * 100
## [1] 2082.9
kbl(summarySE(DryMassExp, measurevar="TotalDryMass", groupvars=c("Aboveground", "Species"),
na.rm = TRUE)[c(2,4),-c(3,5,7)], format = "html",
table.attr = "style='width:75%;'") %>% kable_classic()
| Aboveground | Species | TotalDryMass | se | |
|---|---|---|---|---|
| 2 | Covered | Sapindus saponaria | 0.75867 | 0.16131 |
| 4 | Exposed | Sapindus saponaria | 22.54062 | 12.91268 |
(22.540625 - 0.758667) / 0.758667 * 100
## [1] 2871.1
DryMassExp <- DryMassExp %>%
filter(!is.na(Treatment))
kbl(summarySE(DryMassExp, measurevar="TotalDryMass", groupvars=c("Species", "Treatment"),
na.rm = TRUE)[,-c(5,7)], format = "html",
table.attr = "style='width:75%;'") %>% kable_classic()
| Species | Treatment | N | TotalDryMass | se |
|---|---|---|---|---|
| Cymbopetalum baillonii | No Competition | 3 | 3.17000 | 2.16246 |
| Cymbopetalum baillonii | Belowground Competition | 3 | 4.16000 | 3.11760 |
| Cymbopetalum baillonii | Aboveground Competition | 7 | 0.69571 | 0.14824 |
| Cymbopetalum baillonii | Full Competition | 6 | 1.01000 | 0.33087 |
| Sapindus saponaria | No Competition | 9 | 11.91222 | 5.55222 |
| Sapindus saponaria | Belowground Competition | 7 | 36.20571 | 29.06216 |
| Sapindus saponaria | Aboveground Competition | 8 | 0.61750 | 0.22926 |
| Sapindus saponaria | Full Competition | 7 | 0.92000 | 0.22778 |
Answer: For C. baillonii, the
Only Belowground Competition treatment yielded the highest
dry mass (4.16 ± 3.118), followed by No Competition (31.7 ±
2.162), Full Competition (1.01 ± 0.331) and
Only Aboveground Competition (0.696 ± 0.148).
For S. saponaria, being Belowground Competition
brought the greatest plant mass, but with incredible variation (36.21 ±
29.06), followed by No Competition (11.912 ± 5.552). Lowest
plant mass was a near tie between Full Competition (0.92 ±
0.23) and Aboveground Competition (0.62 ± 0.23).
Did the root to shoot ratio differ between species and treatment?
WetMassExp$Ratio <- WetMassExp$Roots..g. / WetMassExp$Stem..g.
kbl(summarySE(WetMassExp, measurevar="Ratio", groupvars=c("Species", "Treatment"),
na.rm = TRUE)[,-c(5,7)], format = "html",
table.attr = "style='width:75%;'") %>% kable_classic()
| Species | Treatment | N | Ratio | se |
|---|---|---|---|---|
| Cymbopetalum baillonii | No Competition | 3 | 0.57010 | 0.23989 |
| Cymbopetalum baillonii | Belowground Competition | 3 | 0.57595 | 0.25798 |
| Cymbopetalum baillonii | Aboveground Competition | 7 | 0.49621 | 0.16963 |
| Cymbopetalum baillonii | Full Competition | 6 | 0.44570 | 0.13229 |
| Sapindus saponaria | No Competition | 9 | 1.40065 | 0.15759 |
| Sapindus saponaria | Belowground Competition | 7 | 1.13823 | 0.14542 |
| Sapindus saponaria | Aboveground Competition | 8 | 0.84553 | 0.21090 |
| Sapindus saponaria | Full Competition | 7 | 0.80433 | 0.05283 |
Answer: The root-to-shoot ratio in C.
baillonii is highest Only Belowground Competition
(0.576 ± 0.258) and No Competition (0.57 ± 0.239), followed
by Only Aboveground Competition (0.496 ± 0.169), and lastly
Full Competition (0.446 ± 0.132). For S.
saponaria, the root-to-shoot ratio is highest under
No Competition (1.4 ± 0.158), followed by
Only Belowground Competition (1.138 ± 0.145),
Only Aboveground Competition (0.846 ± 0.21), and
Full Competition (0.804 ± 0.053).
Did the dry root to shoot ratio differ between species and treatment?
DryMassExp$Ratio <- DryMassExp$Roots..g. / DryMassExp$Stem..g.
kbl(summarySE(DryMassExp, measurevar="Ratio", groupvars=c("Species", "Treatment"),
na.rm = TRUE)[,-c(5,7)], format = "html",
table.attr = "style='width:75%;'") %>% kable_classic()
| Species | Treatment | N | Ratio | se |
|---|---|---|---|---|
| Cymbopetalum baillonii | No Competition | 3 | 1.08030 | 0.21302 |
| Cymbopetalum baillonii | Belowground Competition | 3 | 0.60492 | 0.20790 |
| Cymbopetalum baillonii | Aboveground Competition | 7 | 0.63790 | 0.13749 |
| Cymbopetalum baillonii | Full Competition | 6 | 0.60153 | 0.14405 |
| Sapindus saponaria | No Competition | 9 | 1.33623 | 0.16288 |
| Sapindus saponaria | Belowground Competition | 7 | 1.11215 | 0.14517 |
| Sapindus saponaria | Aboveground Competition | 8 | 1.12454 | 0.22946 |
| Sapindus saponaria | Full Competition | 7 | 0.86132 | 0.05051 |
Answer: The dry root-to-shoot ratio in C.
baillonii was highest under No Competition (1.08 ±
0.213), followed by Only Aboveground Competition (0.637 ±
0.137), Only Belowground Competition (0.605 ± 0.207), and
No Competition (1.08 ± 0.213).
No Competition (1.336 ± 0.163), followed by
Only Aboveground Competition (1.125 ± 0.229),
Only Belowground Competition (1.112 ± 0.145), and lastly
Full Competition (0.861 ± 0.05).
Does mean leaf mass differ between treatments?
kbl(summarySE(WetMassExp, measurevar="LeafAveMass", groupvars=c("Species", "Treatment"),
na.rm = TRUE)[,-c(5,7)], format = "html",
table.attr = "style='width:75%;'") %>% kable_classic()
| Species | Treatment | N | LeafAveMass | se |
|---|---|---|---|---|
| Cymbopetalum baillonii | No Competition | 3 | 0.23167 | 0.09841 |
| Cymbopetalum baillonii | Belowground Competition | 3 | 0.29767 | 0.21320 |
| Cymbopetalum baillonii | Aboveground Competition | 7 | 0.09686 | 0.01413 |
| Cymbopetalum baillonii | Full Competition | 6 | 0.50883 | 0.39868 |
| Sapindus saponaria | No Competition | 9 | 0.46722 | 0.13258 |
| Sapindus saponaria | Belowground Competition | 7 | 0.71114 | 0.32364 |
| Sapindus saponaria | Aboveground Competition | 8 | 0.07888 | 0.02350 |
| Sapindus saponaria | Full Competition | 7 | 0.10057 | 0.02262 |
Answer: This does appear to be the case. For C.
baillonii, the heaviest leaves were found on the
Full Competition group, with great variance, followed by
the Only Belowground Competition group and
No Competition The Aboveground Competition
group had the lightest leaves. For S. saponaria the heaviest
leaves belonged to the Only Belowground Competition group,
with great variance followed by No Competition. Both
Aboveground Competition and Full Competition
groups had the lightest leaves.
Does mean dry leaf mass differ between treatments?
kbl(summarySE(DryMassExp, measurevar="LeafAveMass", groupvars=c("Species", "Treatment"),
na.rm = TRUE)[,-c(5,7)], format = "html",
table.attr = "style='width:75%;'") %>% kable_classic()
| Species | Treatment | N | LeafAveMass | se |
|---|---|---|---|---|
| Cymbopetalum baillonii | No Competition | 3 | 0.07333 | 0.03930 |
| Cymbopetalum baillonii | Belowground Competition | 3 | 0.11667 | 0.09207 |
| Cymbopetalum baillonii | Aboveground Competition | 7 | 0.03000 | 0.00488 |
| Cymbopetalum baillonii | Full Competition | 6 | 0.03333 | 0.00882 |
| Sapindus saponaria | No Competition | 9 | 0.19000 | 0.05196 |
| Sapindus saponaria | Belowground Competition | 7 | 0.27000 | 0.09061 |
| Sapindus saponaria | Aboveground Competition | 8 | 0.03250 | 0.00840 |
| Sapindus saponaria | Full Competition | 7 | 0.04000 | 0.00900 |
Answer: For C. baillonii, the heaviest
leaves were found on the Only Belowground Competition
treatment (0.117 ± 0.09), followed by No Competition (0.073
± 0.039), Full Competition (0.033 ± 0.009), and lastly
Only Aboveground Competition (0.03 ± 0.005).
For S. saponaria the heaviest leaves belonged to the
Only Belowground Competition treatment (0.27 ± 0.09),
followed by No Competition (0.19 ± 0.519),
Full Competition (0.04 ± 0.009), and lastly
Only Aboveground Competition (0.0325 ± 0.008).
For this, we’ll first need to join our DryMassExp and
WetMassExp dataframes.
DryMassExp$TotalWetMass <- with(WetMassExp, TotalWetMass[match(WetMassExp$TrueID, TrueID)])
DryMassExp$Moisture <- DryMassExp$TotalDryMass / DryMassExp$TotalWetMass * 100
As an alternative, we will recalculate Moisture without the root mass, as the roots were rinsed. We will also do only leaves.
DryMassExp$WetRoots <- with(WetMassExp, Roots..g.[match(WetMassExp$TrueID, TrueID)])
DryMassExp$WetLeaves <- with(WetMassExp, TotalLeafMass[match(WetMassExp$TrueID, TrueID)])
DryMassExp$Moisture_NoRoots <- (DryMassExp$TotalDryMass - DryMassExp$Roots..g.) /
(DryMassExp$TotalWetMass - DryMassExp$WetRoots) * 100
DryMassExp$Moisture_Leaves <- DryMassExp$TotalLeafMass / DryMassExp$WetLeaves * 100
range(DryMassExp$Moisture - DryMassExp$Moisture_NoRoots)
## [1] -4.846 21.775
kbl(summarySE(DryMassExp, measurevar="Moisture", groupvars=c("Species", "Treatment"),
na.rm = TRUE)[,-c(5,7)], format = "html",
table.attr = "style='width:75%;'") %>% kable_classic()
| Species | Treatment | N | Moisture | se |
|---|---|---|---|---|
| Cymbopetalum baillonii | No Competition | 3 | 41.182 | 5.9497 |
| Cymbopetalum baillonii | Belowground Competition | 3 | 43.428 | 4.3983 |
| Cymbopetalum baillonii | Aboveground Competition | 7 | 34.989 | 2.3844 |
| Cymbopetalum baillonii | Full Competition | 6 | 24.926 | 4.2185 |
| Sapindus saponaria | No Competition | 9 | 53.330 | 2.5040 |
| Sapindus saponaria | Belowground Competition | 7 | 58.342 | 3.7171 |
| Sapindus saponaria | Aboveground Competition | 8 | 49.651 | 3.3332 |
| Sapindus saponaria | Full Competition | 7 | 49.023 | 2.8023 |
kbl(summarySE(DryMassExp, measurevar="Moisture_NoRoots", groupvars=c("Species", "Treatment"),
na.rm = TRUE)[,-c(5,7)], format = "html",
table.attr = "style='width:75%;'") %>% kable_classic()
| Species | Treatment | N | Moisture_NoRoots | se |
|---|---|---|---|---|
| Cymbopetalum baillonii | No Competition | 3 | 34.107 | 3.7345 |
| Cymbopetalum baillonii | Belowground Competition | 3 | 41.571 | 4.1753 |
| Cymbopetalum baillonii | Aboveground Competition | 7 | 33.236 | 3.1217 |
| Cymbopetalum baillonii | Full Competition | 6 | 22.725 | 3.9627 |
| Sapindus saponaria | No Competition | 9 | 50.092 | 2.2798 |
| Sapindus saponaria | Belowground Competition | 7 | 56.529 | 4.8130 |
| Sapindus saponaria | Aboveground Competition | 8 | 44.896 | 3.0256 |
| Sapindus saponaria | Full Competition | 7 | 46.151 | 2.4996 |
kbl(summarySE(DryMassExp, measurevar="Moisture_Leaves", groupvars=c("Species", "Treatment"),
na.rm = TRUE)[,-c(5,7)], format = "html",
table.attr = "style='width:75%;'") %>% kable_classic()
| Species | Treatment | N | Moisture_Leaves | se |
|---|---|---|---|---|
| Cymbopetalum baillonii | No Competition | 3 | 29.640 | 2.7795 |
| Cymbopetalum baillonii | Belowground Competition | 3 | 41.387 | 10.7053 |
| Cymbopetalum baillonii | Aboveground Competition | 7 | 31.453 | 4.0788 |
| Cymbopetalum baillonii | Full Competition | 6 | 20.303 | 4.0461 |
| Sapindus saponaria | No Competition | 9 | 41.417 | 2.1253 |
| Sapindus saponaria | Belowground Competition | 7 | 49.186 | 8.0456 |
| Sapindus saponaria | Aboveground Competition | 8 | 40.727 | 4.0108 |
| Sapindus saponaria | Full Competition | 7 | 39.672 | 2.7746 |
Answer: Indeed. For C. baillonii, percent
moisture content was highest in Belowground Competition
(43.43 ± 4.398) and No Competition (41.182 ± 5.949),
followed by Only Aboveground Competition (34.989 ± 2.38),
and lastly Full Competition (24.926 ± 4.219). For S.
saponaria, percent moisture content was highest in
Only Belowground Competition (58.342 ± 3.717), followed by
No Competition (53.33 ± 2.50), and a near tie between
Only Aboveground Competition (49.651 ± 3.333) and
Full Competition (49.023 ± 2.802).
This description is for the table using Moisture not
Moisture_NoRoots or Moisture_Leaves.
Now, lets take a look at the distribution of values for
Count. For this we will create a histogram and kernel
density estimate, which is basically a smoothed out histogram using the
geom_histogram and geom_density arguments of
ggplot. We will also use geom_vline to create
a vertical line where the mean of our species richness is.
ggplot(T7, aes(x=Count)) +
geom_histogram(aes(y=..density..), colour="black", fill="cornsilk") +
geom_density(alpha=.2, fill="#FF6666") + # easier to see distribution
geom_vline(aes(xintercept=mean(Count, na.rm=TRUE)),color="cadetblue4",
linetype="dashed", size=1) + # This adds a blue line where the mean is.
xlab("Number of Germinations") +
ylab("Density") +
theme_bw(base_size = 15)
Interpretation: Our Count data presents a right-skew fueled by some outlier values. The mean is visibly to the right of the peak density of values. This suggests that a normal distribution is unlikely the best way to go for our model.
Still, lets start by making a linear model. We want to find out if there is a difference between species and whether there is an interaction between the experimental treatments. We’ll visualize
options(digits = 6)
lm1 <- lm(data = T7B, Count ~ Species + Aboveground * Belowground)
plot(lm1)
Interpretation: Ok, starting with our residuals vs fitted, we can see that the residuals do not follow the horizontal line all too well. This indicates possible problems with heteroscedasticity and non-normality. The QQ plot shows that we do not have normality of residuals. Our Scale-Location plot does not show residuals with a constant variance, so that confirms heteroscedasticity too.
The Poisson distribution is commonly used when modeling count data, where the variance is equal to the mean. However, if the data show greater variability than expected under a Poisson distribution (i.e. overdispersion), the Negative Binomial distribution can be used as an alternative.
mean(T7B$Count); var(T7B$Count)
## [1] 2.70833
## [1] 3.78546
The variance is larger than the mean, which suggests some
overdispersion. Therefore, a Negative Binomial distribution may be a
good candidate for modeling your data, since it can account for
overdispersion in count data. We can create one using the
glm.nb function of the MASS package.
glm.nb1 <- glm.nb(data = T7B, Count ~ Species * Aboveground * Belowground)
plot(glm.nb1,1); plot(glm.nb1,3)
Interpretation: Our Residuals vs Fitted graph looks better, as does Scale-Location, but I wonder if they could be better.
Let’s perform the Pearson chi-square goodness-of-fit test to determine if the negative binomial is a good fit for our data.
observed_counts <- T7B$Count
expected_counts <- predict(glm.nb1, type = "response")
chisq.test(observed_counts, p = expected_counts / sum(expected_counts))
##
## Chi-squared test for given probabilities
##
## data: observed_counts
## X-squared = 46.76, df = 47, p-value = 0.482
Interpretation: The Negative-Binomial model is a good fit for our germination model (X2 = 46.76, p = 0.482).
Now, lets take a look at the distribution of values for
ProportionDead For this we will create a histogram and
kernel density estimate, which is basically a smoothed out histogram
using the geom_histogram and geom_density
arguments of ggplot. We will also use
geom_vline to create a vertical line where the mean of our
species richness is.
ggplot(T12B, aes(x=ProportionDead)) +
geom_histogram(aes(y=..density..), colour="black", fill="cornsilk") +
geom_density(alpha=.2, fill="#FF6666") + # easier to see distribution
geom_vline(aes(xintercept=mean(ProportionDead, na.rm=TRUE)),color="cadetblue4",
linetype="dashed", size=1) + # This adds a blue line where the mean is.
xlab("Proportion of Deaths") +
ylab("Density") +
theme_bw(base_size = 15)
Interpretation: Well, our plot looks a little bimodal! High zeroes and high ones, both about equal, so the density graph is flat. Not sure what kind of model structure best suits this type of distribution.
Let’s start with a linear model.
lm2 <- lm(data = T12B, ProportionDead ~ Species * Aboveground * Belowground)
plot(lm2)
Interpretation: Residuals vs. Fitted looks great, and the QQ Plot does not look too shabby either. Scale-Location looks wack though, yo. 85 and 91 appear to be the consistently high value blocks, but again, no reason to go outlier hunting. Let’s try our test for Poisson and Negative Binomial.
mean(T12B$ProportionDead); var(T12B$ProportionDead)
## [1] 0.515455
## [1] 0.151667
glm.nb2 <- glm.nb(data = T12B, ProportionDead ~ Species * Aboveground * Belowground)
glm.nb3 <- glm.nb(data = T12B, ProportionDead ~ Species + Aboveground * Belowground)
glm.nb4 <- glm.nb(data = T12B, ProportionDead ~ Species + Aboveground + Belowground)
plot(glm.nb2,1); plot(glm.nb2,3)
observed_counts <- T12B$Count
expected_counts <- predict(glm.nb2, type = "response")
chisq.test(observed_counts, p = expected_counts / sum(expected_counts))
##
## Chi-squared test for given probabilities
##
## data: observed_counts
## X-squared = 32.65, df = 43, p-value = 0.874
Interpretation: Poisson is unlikely to be a good fit given the difference between mean and variance, but Negative Binomial does pass the chi-square test though (X2 = 30.665, p = 0.8808), and the plots looks good enough.
lm3 <- lm(data = FernExp, Weeks.Alive ~ Species * Aboveground * Belowground)
plot(lm3)
Interpretation: Our Residuals vs Fittted plot looks great, but the QQ plot is very S-shaped. Let’s check if Poisson works for our data.
mean(FernExp$Weeks.Alive, na.rm = TRUE)
## [1] 44.2205
var(FernExp$Weeks.Alive, na.rm = TRUE)
## [1] 1451.11
Interpretation: We have overdispersion, so Poisson probably would not be a good idea. Let’s build our Negative-Binomial model once again and check assumptions.
FernExp$Belowground <- relevel(FernExp$Belowground, ref = "Unprotected")
glm.nb5 <- glm.nb(data = FernExp, Weeks.Alive ~ Species + Aboveground * Belowground)
glm.nb6 <- glm.nb(data = FernExp, Weeks.Alive ~ Species + Aboveground * Belowground)
glm.nb7 <- glm.nb(data = FernExp, Weeks.Alive ~ Species + Aboveground + Belowground)
plot(glm.nb5,1); plot(glm.nb5,3)
observed_counts <- na.omit(FernExp$Weeks.Alive)
expected_counts <- predict(glm.nb5, type = "response")
chisq.test(observed_counts, p = expected_counts / sum(expected_counts))
##
## Chi-squared test for given probabilities
##
## data: observed_counts
## X-squared = 4042, df = 126, p-value <2e-16
Interpretation: According to the Chi-Squared Test, the negative binomial is not a good fit either, despite the diagnostic graphs looking fairly well. Let’s try Gamma.
glm.gamma.1 <- glm(data = FernExp, Weeks.Alive ~ Species * Aboveground * Belowground,
family = Gamma(link = "log"))
glm.gamma.2 <- glm(data = FernExp, Weeks.Alive ~ Species + Aboveground * Belowground,
family = Gamma(link = "log"))
glm.gamma.3 <- glm(data = FernExp, Weeks.Alive ~ Species + Aboveground + Belowground,
family = Gamma(link = "log"))
plot(glm.gamma.1)
Interpretation: Gamma looks a little better in my opinion. We’ll stick with Gamma.
ggplot(FernExp.G2, aes(x=MaxLeaves)) +
geom_histogram(aes(y=..density..), colour="black", fill="cornsilk") +
geom_density(alpha=.2, fill="#FF6666") + # easier to see distribution
geom_vline(aes(xintercept=mean(MaxLeaves, na.rm=TRUE)),color="cadetblue4",
linetype="dashed", size=1) + # This adds a blue line where the mean is.
xlab("Max Leaf Count") +
ylab("Density") +
theme_bw(base_size = 15)
Interpretation: Welp, another right-skewed distribution. Who said negative-binomial forever? Let’s of course first check for Poisson. We might need something that is zero-inflated.
mean(FernExp.G2$MaxLeaves)
## [1] 5.8
var(FernExp.G2$MaxLeaves)
## [1] 44.2258
Interpretation: Not even close – this variable is even more overdispersed. Let’s create the negative binomial model. Different species are expected to have different leaf counts, so let’s create separate models for each species.
glm.nb8.CYBA <- glm.nb(data = FernExp.G.CYBA, MaxLeaves ~ Aboveground * Belowground)
glm.nb8.SASA <- glm.nb(data = FernExp.G.SASA, MaxLeaves ~ Aboveground * Belowground)
glm.nb8 <- glm.nb(data = FernExp.G2, MaxLeaves ~ Species * Aboveground * Belowground)
observed_counts <- FernExp.G.CYBA$MaxLeaves
expected_counts <- predict(glm.nb8.CYBA, type = "response")
chisq.test(observed_counts, p = expected_counts / sum(expected_counts))
##
## Chi-squared test for given probabilities
##
## data: observed_counts
## X-squared = 539.6, df = 77, p-value <2e-16
observed_counts <- FernExp.G.SASA$MaxLeaves
expected_counts <- predict(glm.nb8.SASA, type = "response")
chisq.test(observed_counts, p = expected_counts / sum(expected_counts))
##
## Chi-squared test for given probabilities
##
## data: observed_counts
## X-squared = 250.5, df = 46, p-value <2e-16
observed_counts <- FernExp.G2$MaxLeaves
expected_counts <- predict(glm.nb8, type = "response")
chisq.test(observed_counts, p = expected_counts / sum(expected_counts))
##
## Chi-squared test for given probabilities
##
## data: observed_counts
## X-squared = 790.1, df = 124, p-value <2e-16
Interpretation: Ah, well, the negative binomial Chi-Square test check is giving very small p-values. This might have to do with that zero-inflation. Let’s look at some basic stats.
summary(FernExp.G2$MaxLeaves)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0 2.0 3.0 5.8 8.0 36.0
We have some zeroes, but mostly a very right-skewed distribution. Let’s check how many zeroes we have.
nrow(droplevels(subset(FernExp.G2, MaxLeaves == "0"))) / nrow(FernExp.G2) * 100
## [1] 21.6
Ah, a fifth are zeroes. Perhaps we should use a zero-inflated model.
Let’s try a Zero-Inflated Negative Binomial distribution (ZINB),
achievable using the zeroinfl function in the
pscl package. The ZINB distribution allows for
overdispersion (i.e., variance is greater than the mean) in count data
and models both the probability of extra zeros and the variance of the
non-zero counts. We can then use the lrtest function to verify that the
model fits better than a regular negative binomial.
zinb_CYBA <- zeroinfl(MaxLeaves ~ Aboveground * Belowground | 1, data = FernExp.G.CYBA, dist = "negbin")
zinb_SASA <- zeroinfl(MaxLeaves ~ Aboveground * Belowground | 1, data = FernExp.G.SASA, dist = "negbin")
zinb_2 <- zeroinfl(MaxLeaves ~ Species * Aboveground * Belowground | 1, data = FernExp.G2, dist = "negbin")
lrtest(glm.nb8.CYBA, zinb_CYBA)
## Likelihood ratio test
##
## Model 1: MaxLeaves ~ Aboveground * Belowground
## Model 2: MaxLeaves ~ Aboveground * Belowground | 1
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 5 -199.8
## 2 6 -199.1 1 1.282 0.258
lrtest(glm.nb8.SASA, zinb_SASA)
## Likelihood ratio test
##
## Model 1: MaxLeaves ~ Aboveground * Belowground
## Model 2: MaxLeaves ~ Aboveground * Belowground | 1
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 5 -144.3
## 2 6 -139.9 1 8.791 0.00303 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lrtest(glm.nb8, zinb_2)
## Likelihood ratio test
##
## Model 1: MaxLeaves ~ Species * Aboveground * Belowground
## Model 2: MaxLeaves ~ Species * Aboveground * Belowground | 1
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 9 -347.1
## 2 10 -343.1 1 8.085 0.00446 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Interpretation: Indeed, the Likelihood Ratio Tests show that the ZINB models are superior. Let’s make variants with random effect of plot.
zinb_SASA_GLMMTMB <- glmmTMB(MaxLeaves ~ Aboveground * Belowground + (1 | Plot),
ziformula = ~ 1,
family = nbinom2(),
data = FernExp.G.SASA)
zinb_CYBA_GLMMTMB <- glmmTMB(MaxLeaves ~ Aboveground * Belowground + (1 | Plot),
ziformula = ~ 1,
family = nbinom2(),
data = FernExp.G.CYBA)
ggplot(WetMassExp, aes(x=Height..cm.)) +
geom_histogram(aes(y=..density..), colour="black", fill="cornsilk") +
geom_density(alpha=.2, fill="#FF6666") + # easier to see distribution
geom_vline(aes(xintercept=mean(Height..cm., na.rm=TRUE)),color="cadetblue4",
linetype="dashed", size=1) + # This adds a blue line where the mean is.
xlab("Seedling Height on Harvest") +
ylab("Density") +
theme_bw(base_size = 15)
Interpretation: We have a right-skewed distribution of positive numbers, so lets go with Gamma again.
glmHeight <- glm(data = WetMassExp, Height..cm. ~ Species * Aboveground * Belowground,
family = Gamma(link = "log"))
plot(glmHeight)
Interpretation: Our Gamma looks good! The outlier we see is expected. Just an unusually tall plant.
ggplot(WetMassExp, aes(x=Diameter.Base..cm.)) +
geom_histogram(aes(y=..density..), colour="black", fill="cornsilk") +
geom_density(alpha=.2, fill="#FF6666") + # easier to see distribution
geom_vline(aes(xintercept=mean(LeafAveMass, na.rm=TRUE)),color="cadetblue4",
linetype="dashed", size=1) + # This adds a blue line where the mean is.
xlab("Base Diameter (cm)") +
ylab("Density") +
theme_bw(base_size = 15)
Very right-skewed with a hell of an outlier, likely in that S. saponaria Exposed and Unprotected category. Let’s look at this average again without the outlier.
WetMassExp.C <- droplevels(subset(WetMassExp, WetMassExp$Diameter.Base..cm. < 9))
kbl(summarySE(WetMassExp.C, measurevar="Diameter.Base..cm.", groupvars=c("Species", "Treatment"),
na.rm = TRUE)[,-c(5,7)], format = "html",
table.attr = "style='width:75%;'") %>% kable_classic()
| Species | Treatment | N | Diameter.Base..cm. | se |
|---|---|---|---|---|
| Cymbopetalum baillonii | No Competition | 3 | 0.433333 | 0.088192 |
| Cymbopetalum baillonii | Belowground Competition | 3 | 0.666667 | 0.272845 |
| Cymbopetalum baillonii | Aboveground Competition | 7 | 0.357143 | 0.020203 |
| Cymbopetalum baillonii | Full Competition | 6 | 0.350000 | 0.042817 |
| Sapindus saponaria | No Competition | 9 | 0.722222 | 0.105116 |
| Sapindus saponaria | Belowground Competition | 6 | 0.950000 | 0.416133 |
| Sapindus saponaria | Aboveground Competition | 8 | 0.287500 | 0.029505 |
| Sapindus saponaria | Full Competition | 7 | 0.328571 | 0.035952 |
ggplot(WetMassExp.C, aes(x=Diameter.Base..cm.)) +
geom_histogram(aes(y=..density..), colour="black", fill="cornsilk") +
geom_density(alpha=.2, fill="#FF6666") + # easier to see distribution
geom_vline(aes(xintercept=mean(LeafAveMass, na.rm=TRUE)),color="cadetblue4",
linetype="dashed", size=1) + # This adds a blue line where the mean is.
xlab("Base Diameter (cm)") +
ylab("Density") +
theme_bw(base_size = 15)
Interpretation: Removing the outlier lowers the mean and se considerably. There is still a right-skew, but it is not as extreme. Should we be removing this mega outlier?
glmG.Diam <- glm(data = WetMassExp, Diameter.Base..cm. ~ Species * Aboveground * Belowground,
family = Gamma(link = "log"))
glmG.DiamEX <- glm(data = WetMassExp.C, Diameter.Base..cm. ~ Species * Aboveground * Belowground,
family = Gamma(link = "log"))
plot(glmG.Diam)
Interpretation: Gamma looks fine, but that extreme outlier is a bit of a concern. Let’s run both versions and check the differences.
First we’ll want to check the distribution of our data, but what is immediately apparent, and something I remember from the harvest, is that we have one HUGE plant of S. saponaria, we’re going to want to visualize what the distribution looks like with and without this giant.
WetMassExpB <- droplevels(subset(WetMassExp, WetMassExp$TotalWetMass < 200))
ggplot(WetMassExp, aes(x=TotalWetMass)) +
geom_histogram(aes(y=..density..), colour="black", fill="cornsilk") +
geom_density(alpha=.2, fill="#FF6666") + # easier to see distribution
geom_vline(aes(xintercept=mean(TotalWetMass, na.rm=TRUE)),color="cadetblue4",
linetype="dashed", size=1) + # This adds a blue line where the mean is.
xlab("Total Fresh Mass (g)") +
ylab("Density") +
theme_bw(base_size = 15)
ggplot(WetMassExpB, aes(x=TotalWetMass)) +
geom_histogram(aes(y=..density..), colour="black", fill="cornsilk") +
geom_density(alpha=.2, fill="#FF6666") + # easier to see distribution
geom_vline(aes(xintercept=mean(TotalWetMass, na.rm=TRUE)),color="cadetblue4",
linetype="dashed", size=1) + # This adds a blue line where the mean is.
xlab("Total Fresh Mass (g)") +
ylab("Density") +
theme_bw(base_size = 15)
Ok! So, regardless of whether we include the outlier, we have a right-skewed distribution. A normal model is unlikely to work, but perhaps we can use a Gamma, considering we do not have zeroes. Let’s check the normal model.
lm4 <- lm(data = WetMassExpB, TotalWetMass ~ Species * Aboveground * Belowground)
plot(lm4)
Interpretation: Yeah, the Residuals vs Fitted looks barely okay, the QQ-Plot has a jump in values, which is what we expect to see but still not ideal for a normal model, let’s just go ahead and try Gamma.
glmFresh <- glm(data = WetMassExp, TotalWetMass ~ Species * Aboveground * Belowground,
family = Gamma(link = "log"))
glmFresh2 <- glm(data = WetMassExp, TotalWetMass ~ Species + Aboveground * Belowground,
family = Gamma(link = "log"))
glmFresh3 <- glm(data = WetMassExp, TotalWetMass ~ Species + Aboveground + Belowground,
family = Gamma(link = "log"))
plot(glmFresh)
Interpretation: The Residual vs Fitted and QQ Plots look better! We still have values curving away, but given the treatment effects, this is to be expected. Also, despite including our super high value, it does not pass the Cook’s distance.
First we’ll want to check the distribution of our data, but what is immediately apparent, and something I remember from the harvest, is that we have one HUGE plant of S. saponaria, we’re going to want to visualize what the distribution looks like with and without this giant.
DryMassExpB <- droplevels(subset(DryMassExp, DryMassExp$TotalDryMass < 200))
ggplot(DryMassExp, aes(x=TotalDryMass)) +
geom_histogram(aes(y=..density..), colour="black", fill="cornsilk") +
geom_density(alpha=.2, fill="#FF6666") + # easier to see distribution
geom_vline(aes(xintercept=mean(TotalDryMass, na.rm=TRUE)),color="cadetblue4",
linetype="dashed", size=1) + # This adds a blue line where the mean is.
xlab("Total Fresh Mass (g)") +
ylab("Density") +
theme_bw(base_size = 15)
ggplot(DryMassExpB, aes(x=TotalDryMass)) +
geom_histogram(aes(y=..density..), colour="black", fill="cornsilk") +
geom_density(alpha=.2, fill="#FF6666") + # easier to see distribution
geom_vline(aes(xintercept=mean(TotalDryMass, na.rm=TRUE)),color="cadetblue4",
linetype="dashed", size=1) + # This adds a blue line where the mean is.
xlab("Total Fresh Mass (g)") +
ylab("Density") +
theme_bw(base_size = 15)
Ok! So, regardless of whether we include the outlier, we have a right-skewed distribution. A normal model is unlikely to work, but perhaps we can use a Gamma, considering we do not have zeroes. Let’s check the normal model.
lm4 <- lm(data = DryMassExpB, TotalDryMass ~ Species * Aboveground * Belowground)
plot(lm4)
Interpretation: Yeah, the Residuals vs Fitted looks barely okay, the QQ-Plot is not great and shows the same jump in values, which is what we expect to see but still not ideal for a normal model, let’s just go ahead and try Gamma.
DryMassExp$Belowground <- relevel(DryMassExp$Belowground, ref = "Unprotected")
glmDry <- glm(data = DryMassExp, TotalDryMass ~ Species * Aboveground * Belowground,
family = Gamma(link = "log"))
plot(glmDry)
Interpretation: The Residual vs Fitted and QQ Plots look better! We still have values curving away, but given the treatment effects, this is to be expected. Also, despite including our super high value, it does not pass the Cook’s distance.
ggplot(WetMassExp, aes(x=Ratio)) +
geom_histogram(aes(y=..density..), colour="black", fill="cornsilk") +
geom_density(alpha=.2, fill="#FF6666") + # easier to see distribution
geom_vline(aes(xintercept=mean(Ratio, na.rm=TRUE)),color="cadetblue4",
linetype="dashed", size=1) + # This adds a blue line where the mean is.
xlab("Root-to-Shoot Ratio") +
ylab("Density") +
theme_bw(base_size = 15)
Interpretation: There’s a chance that we can use the normal gaussian on this.
lm5Wet <- lm(data = WetMassExp, Ratio ~ Species * Aboveground * Belowground)
plot(lm5Wet)
Interpretation: Yes, it does seem like the gaussian distribution is a good enough fit.
ggplot(DryMassExp, aes(x=Ratio)) +
geom_histogram(aes(y=..density..), colour="black", fill="cornsilk") +
geom_density(alpha=.2, fill="#FF6666") + # easier to see distribution
geom_vline(aes(xintercept=mean(Ratio, na.rm=TRUE)),color="cadetblue4",
linetype="dashed", size=1) + # This adds a blue line where the mean is.
xlab("Dry Root-to-Shoot Ratio") +
ylab("Density") +
theme_bw(base_size = 15)
Interpretation: Looks normal-ish!
lm5Dry <- lm(data = DryMassExp, Ratio ~ Species * Aboveground * Belowground)
plot(lm5Dry)
Interpretation: Graphs looking a little ok. The QQ-Plot shows treatment effects. The Scale-Location Plot is the most odd, but lets try modeling with a normal model anyway.
ggplot(WetMassExp, aes(x=LeafAveMass)) +
geom_histogram(aes(y=..density..), colour="black", fill="cornsilk") +
geom_density(alpha=.2, fill="#FF6666") + # easier to see distribution
geom_vline(aes(xintercept=mean(LeafAveMass, na.rm=TRUE)),color="cadetblue4",
linetype="dashed", size=1) + # This adds a blue line where the mean is.
xlab("Mean Leaf Mass (g)") +
ylab("Density") +
theme_bw(base_size = 15)
Very right-skewed. We might be able to use Gamma for this one since we have no zeroes.
WetMassExp.SASA <- droplevels(subset(WetMassExp, WetMassExp$Species == "Sapindus saponaria"))
WetMassExp.CYBA <- droplevels(subset(WetMassExp, WetMassExp$Species == "Cymbopetalum baillonii"))
lm8.CYBA <- lm(data = WetMassExp.CYBA, LeafAveMass ~ Aboveground * Belowground)
lm8.SASA <- lm(data = WetMassExp.SASA, LeafAveMass ~Aboveground * Belowground)
plot(lm8.CYBA)
plot(lm8.SASA)
Interpretation: The CYBA and SASA models look as expected. The outlier is expected. Linear model not great. Let’s check to see if Poisson could work.
mean(WetMassExp.CYBA$LeafAveMass); var(WetMassExp.CYBA$LeafAveMass)
## [1] 0.279947
## [1] 0.314696
mean(WetMassExp.SASA$LeafAveMass); var(WetMassExp.SASA$LeafAveMass)
## [1] 0.33929
## [1] 0.259131
mean(WetMassExp$LeafAveMass); var(WetMassExp$LeafAveMass)
## [1] 0.31674
## [1] 0.2751
Some underdispersion in both cases.
glm.nb.Leaf.CYBA <- glm.nb(data = WetMassExp.CYBA, LeafAveMass ~ Aboveground * Belowground)
glm.nb.Leaf.SASA <- glm.nb(data = WetMassExp.SASA, LeafAveMass ~ Aboveground * Belowground)
glm.nb.Leaf <- glm.nb(data = WetMassExp, LeafAveMass ~ Species * Aboveground * Belowground)
observed_counts <- WetMassExp.CYBA$LeafAveMass
expected_counts <- predict(glm.nb.Leaf.CYBA, type = "response")
chisq.test(observed_counts, p = expected_counts / sum(expected_counts))
##
## Chi-squared test for given probabilities
##
## data: observed_counts
## X-squared = 10.62, df = 18, p-value = 0.91
observed_counts <- WetMassExp.SASA$LeafAveMass
expected_counts <- predict(glm.nb.Leaf.SASA, type = "response")
chisq.test(observed_counts, p = expected_counts / sum(expected_counts))
##
## Chi-squared test for given probabilities
##
## data: observed_counts
## X-squared = 9.501, df = 30, p-value = 1
observed_counts <- WetMassExp$LeafAveMass
expected_counts <- predict(glm.nb.Leaf, type = "response")
chisq.test(observed_counts, p = expected_counts / sum(expected_counts))
##
## Chi-squared test for given probabilities
##
## data: observed_counts
## X-squared = 20.13, df = 49, p-value = 1
Interpretation: Looks like the negative binomial distribution will work nicely for our CYBA and SASA models.
ggplot(DryMassExp, aes(x=LeafAveMass)) +
geom_histogram(aes(y=..density..), colour="black", fill="cornsilk") +
geom_density(alpha=.2, fill="#FF6666") + # easier to see distribution
geom_vline(aes(xintercept=mean(LeafAveMass, na.rm=TRUE)),color="cadetblue4",
linetype="dashed", size=1) + # This adds a blue line where the mean is.
xlab("Mean Leaf Mass (g)") +
ylab("Density") +
theme_bw(base_size = 15)
Very right-skewed. We might be able to use Gamma for this one since we have no zeroes.
DryMassExp$Moisture2 <- DryMassExp$Moisture/100
DryMassExp$Moisture_NoRoots2 <- DryMassExp$Moisture_NoRoots/100
DryMassExp$Moisture_Leaves2 <- DryMassExp$Moisture_Leaves/100
DryMassExp.SASA <- droplevels(subset(DryMassExp, DryMassExp$Species == "Sapindus saponaria"))
DryMassExp.CYBA <- droplevels(subset(DryMassExp, DryMassExp$Species == "Cymbopetalum baillonii"))
lm8.CYBA <- lm(data = DryMassExp.CYBA, LeafAveMass ~ Aboveground * Belowground)
lm8.SASA <- lm(data = DryMassExp.SASA, LeafAveMass ~Aboveground * Belowground)
plot(lm8.CYBA)
plot(lm8.SASA)
Interpretation: The CYBA and SASA models look as expected. The outlier is expected. Linear model not great. Let’s check to see if Poisson could work.
mean(WetMassExp.CYBA$LeafAveMass); var(WetMassExp.CYBA$LeafAveMass)
## [1] 0.279947
## [1] 0.314696
mean(WetMassExp.SASA$LeafAveMass); var(WetMassExp.SASA$LeafAveMass)
## [1] 0.33929
## [1] 0.259131
Some underdispersion in both cases.
glm.nb.Leaf.CYBA <- glm.nb(data = DryMassExp.CYBA, LeafAveMass ~ Aboveground * Belowground)
glm.nb.Leaf.SASA <- glm.nb(data = DryMassExp.SASA, LeafAveMass ~ Aboveground * Belowground)
glm.nb.Leaf <- glm.nb(data = DryMassExp, LeafAveMass ~ Species * Aboveground * Belowground)
observed_counts <- DryMassExp.CYBA$LeafAveMass
expected_counts <- predict(glm.nb.Leaf.CYBA, type = "response")
chisq.test(observed_counts, p = expected_counts / sum(expected_counts))
##
## Chi-squared test for given probabilities
##
## data: observed_counts
## X-squared = 0.6657, df = 18, p-value = 1
observed_counts <- DryMassExp.SASA$LeafAveMass
expected_counts <- predict(glm.nb.Leaf.SASA, type = "response")
chisq.test(observed_counts, p = expected_counts / sum(expected_counts))
##
## Chi-squared test for given probabilities
##
## data: observed_counts
## X-squared = 2.507, df = 30, p-value = 1
observed_counts <- DryMassExp$LeafAveMass
expected_counts <- predict(glm.nb.Leaf, type = "response")
chisq.test(observed_counts, p = expected_counts / sum(expected_counts))
##
## Chi-squared test for given probabilities
##
## data: observed_counts
## X-squared = 3.172, df = 49, p-value = 1
Interpretation: Looks like the negative binomial distribution will work nicely for SASA and CYBA. However, looking at the confidence intervals further along, the intervals for CYBA are unusually large, so we will use a Gamma distribution.
DryMassExp.CYBA.B <- glm(data = DryMassExp.CYBA, LeafAveMass ~ Aboveground * Belowground,
family = Gamma(link = "log"))
plot(DryMassExp.CYBA.B)
DryMassExp.SASA.B <- glm(data = DryMassExp.SASA, LeafAveMass ~ Aboveground * Belowground,
family = Gamma(link = "log"))
plot(DryMassExp.SASA.B)
DryMassExp.B <- glm(data = DryMassExp, LeafAveMass ~ Species * Aboveground * Belowground,
family = Gamma(link = "log"))
plot(DryMassExp.B)
Yeah, let’s stick with negative binomial.
ggplot(DryMassExp, aes(x=Moisture)) +
geom_histogram(aes(y=..density..), colour="black", fill="cornsilk") +
geom_density(alpha=.2, fill="#FF6666") + # easier to see distribution
geom_vline(aes(xintercept=mean(Moisture, na.rm=TRUE)),color="cadetblue4",
linetype="dashed", size=1) + # This adds a blue line where the mean is.
xlab("Moisture Content (%)") +
ylab("Density") +
theme_bw(base_size = 15)
ggplot(DryMassExp, aes(x=Moisture_NoRoots)) +
geom_histogram(aes(y=..density..), colour="black", fill="cornsilk") +
geom_density(alpha=.2, fill="#FF6666") + # easier to see distribution
geom_vline(aes(xintercept=mean(Moisture_NoRoots, na.rm=TRUE)),color="cadetblue4",
linetype="dashed", size=1) + # This adds a blue line where the mean is.
xlab("Moisture Content (%)") +
ylab("Density") +
theme_bw(base_size = 15)
ggplot(DryMassExp, aes(x=Moisture_Leaves)) +
geom_histogram(aes(y=..density..), colour="black", fill="cornsilk") +
geom_density(alpha=.2, fill="#FF6666") +
geom_vline(aes(xintercept=mean(Moisture_Leaves, na.rm=TRUE)),color="cadetblue4",
linetype="dashed", size=1) +
xlab("Moisture Content (%)") +
ylab("Density") +
theme_bw(base_size = 15)
Interpretation: Hmm. Well, we should be able to use
a normal gaussian for this, but the confidence intervals we get later
are huge, which suggests this is not a good idea. Lets use beta
regression instead using the betareg function. Beta
regression is specifically designed for modeling continuous percentages
and can handle data bounded between 0 and 1. It models the response
variable using a beta distribution and employs a suitable link function.
The betareg package in R provides functions to fit beta regression
models. If not, we can always just use the normal gaussian, as this does
look pretty normal to me.
betaM.1 <- betareg(data = DryMassExp, Moisture2 ~ Species * Aboveground * Belowground)
betaM.2 <- betareg(data = DryMassExp, Moisture2 ~ Species + Aboveground * Belowground)
betaM.3 <- betareg(data = DryMassExp, Moisture2 ~ Species + Treatment)
betaM.4 <- betareg(data = DryMassExp, Moisture2 ~ Species)
plot(betaM.1)
Interpretation: Indeed, our diagnostic graphs look fairly good considering what we expect from our dataset. It looks like we still have a higher than usual value there, perhaps the huge S. saponaria did not dry completely?
First let’s make a copy of our data without the shadehouse treatment.
Fern2 <- droplevels(subset(Fern, Fern$Plot != "Shadehouse"))
levels(Fern2$Treatment) <- list("No Competition" = "Exposed + Protected",
"Only Belowground Competition" = "Exposed + Unprotected",
"Only Aboveground Competition" = "Covered + Protected",
"Full Competition" = "Covered + Unprotected")
Kaplan-Meier analysis, also known as the Kaplan-Meier method or the product-limit method, is a non-parametric statistical technique used to estimate the survival probability over time in medical research and other fields where time-to-event data is analyzed.
The first thing to do is to use Surv function to build
the standard survival object. The variable Weeks.Alive
records survival time; Status indicates whether the plant’s
death was observed (status = 1) or that survival time was censored
(status = 0). The survfit function is used to produce the
Kaplan-Meier estimates of the probability of survival over time.
Fern2 <- Fern2 %>%
filter(!(ID %in% exclude_ids))
Fern2.SASA <- droplevels(subset(Fern2, Fern2$Species == "Sapindus saponaria"))
Fern2.CYBA <- droplevels(subset(Fern2, Fern2$Species == "Cymbopetalum baillonii"))
fi_trt_fit_SASA <- survfit(Surv(Weeks.Alive, Status) ~ Treatment, data=Fern2.SASA)
summary(as.factor(Fern2.SASA$Status))
## 0 1 NA's
## 31 15 189
summary(Fern2.SASA$Weeks.Alive)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 1.0 25.0 52.0 54.4 90.0 97.0 190
summary(Fern2.SASA$Treatment)
## No Competition Only Belowground Competition
## 55 60
## Only Aboveground Competition Full Competition
## 60 60
Kaplan_TrT_SASA <- autoplot(fi_trt_fit_SASA) +
ylab("Probability of Survival") +
xlab("Weeks Since Germination") +
geom_label(aes(x = 1, y = 0,
label = paste0("p = 0.7")),
vjust = -9.6, hjust = -5.825,
size = 5, fill = "white", color = "black") +
scale_color_manual(values = c("black", "black", "black", "black"),
guide = guide_legend(direction = "horizontal",
title = NULL,
nrow = 2, ncol = 2)) +
scale_fill_manual(values = c("hotpink", "#1F78B4", "limegreen", "#FFD700"),
guide = guide_legend(direction = "horizontal",
title = NULL,
nrow = 2, ncol = 2)) +
theme_classic(base_size = 15) +
theme(legend.position = "bottom",
legend.background = element_rect(color = "black")) +
guides(color = FALSE)
Kaplan_TrT_SASA
survdiff(Surv(Weeks.Alive, Status) ~ Treatment, data = Fern2.SASA)
## Call:
## survdiff(formula = Surv(Weeks.Alive, Status) ~ Treatment, data = Fern2.SASA)
##
## n=45, 190 observations deleted due to missingness.
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## Treatment=No Competition 13 6 4.30 0.6767 0.9623
## Treatment=Only Belowground Competition 10 3 3.53 0.0806 0.1069
## Treatment=Only Aboveground Competition 10 2 3.62 0.7216 0.9662
## Treatment=Full Competition 12 4 3.56 0.0555 0.0738
##
## Chisq= 1.6 on 3 degrees of freedom, p= 0.7
png("~/Downloads/KaplanTreatment_SASA.png", width = 165, height = 120, units = 'mm', res = 300)
Kaplan_TrT_SASA
dev.off()
## quartz_off_screen
## 2
fi_trt_fit_SASA_AG <- survfit(Surv(Weeks.Alive, Status) ~ Aboveground, data=Fern2.SASA)
Kaplan_TrT_SASA_AG <- autoplot(fi_trt_fit_SASA_AG) +
ylab("Probability of Survival") +
xlab("Weeks Since Germination") +
geom_label(aes(x = 1, y = 0,
label = paste0("p = 0.5")),
vjust = -9.6, hjust = -5.825,
size = 5, fill = "white", color = "black") +
scale_color_manual(values = c("black", "black", "black", "black"),
guide = guide_legend(direction = "horizontal",
title = NULL,
nrow = 2, ncol = 2)) +
scale_fill_manual(values = c("hotpink", "#1F78B4", "limegreen", "#FFD700"),
guide = guide_legend(direction = "horizontal",
title = NULL,
nrow = 2, ncol = 2)) +
theme_classic(base_size = 15) +
theme(legend.position = "bottom",
legend.background = element_rect(color = "black")) +
guides(color = FALSE)
Kaplan_TrT_SASA_AG
survdiff(Surv(Weeks.Alive, Status) ~ Aboveground, data = Fern2.SASA)
## Call:
## survdiff(formula = Surv(Weeks.Alive, Status) ~ Aboveground, data = Fern2.SASA)
##
## n=45, 190 observations deleted due to missingness.
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## Aboveground=Covered 22 6 7.17 0.191 0.372
## Aboveground=Exposed 23 9 7.83 0.175 0.372
##
## Chisq= 0.4 on 1 degrees of freedom, p= 0.5
png("~/Downloads/KaplanAG_SASA.png", width = 165, height = 120, units = 'mm', res = 300)
Kaplan_TrT_SASA_AG
dev.off()
## quartz_off_screen
## 2
fi_trt_fit_SASA_BG <- survfit(Surv(Weeks.Alive, Status) ~ Belowground, data=Fern2.SASA)
Kaplan_TrT_SASA_BG <- autoplot(fi_trt_fit_SASA_BG) +
ylab("Probability of Survival") +
xlab("Weeks Since Germination") +
geom_label(aes(x = 1, y = 0,
label = paste0("p = 0.5")),
vjust = -9.6, hjust = -5.825,
size = 5, fill = "white", color = "black") +
scale_color_manual(values = c("black", "black", "black", "black"),
guide = guide_legend(direction = "horizontal",
title = NULL,
nrow = 2, ncol = 2)) +
scale_fill_manual(values = c("hotpink", "#1F78B4", "limegreen", "#FFD700"),
guide = guide_legend(direction = "horizontal",
title = NULL,
nrow = 2, ncol = 2)) +
theme_classic(base_size = 15) +
theme(legend.position = "bottom",
legend.background = element_rect(color = "black")) +
guides(color = FALSE)
Kaplan_TrT_SASA_BG
survdiff(Surv(Weeks.Alive, Status) ~ Belowground, data = Fern2.SASA)
## Call:
## survdiff(formula = Surv(Weeks.Alive, Status) ~ Belowground, data = Fern2.SASA)
##
## n=45, 190 observations deleted due to missingness.
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## Belowground=Protected 23 8 7.91 0.00102 0.00218
## Belowground=Unprotected 22 7 7.09 0.00113 0.00218
##
## Chisq= 0 on 1 degrees of freedom, p= 1
png("~/Downloads/KaplanBG_SASA.png", width = 165, height = 120, units = 'mm', res = 300)
Kaplan_TrT_SASA_BG
dev.off()
## quartz_off_screen
## 2
fi_trt_fit_CYBA <- survfit(Surv(Weeks.Alive, Status) ~ Treatment, data=Fern2.CYBA)
Kaplan_TrT_CYBA <- autoplot(fi_trt_fit_CYBA) +
ylab("Probability of Survival") +
xlab("Weeks Since Germination") +
geom_label(aes(x = 1, y = 0,
label = paste0("p = 0.5")),
vjust = -9.6, hjust = -5.825,
size = 5, fill = "white", color = "black") +
scale_color_manual(values = c("black", "black", "black", "black"),
guide = guide_legend(direction = "horizontal",
title = NULL,
nrow = 2, ncol = 2)) +
scale_fill_manual(values = c("hotpink", "#1F78B4", "limegreen", "#FFD700"),
guide = guide_legend(direction = "horizontal",
title = NULL,
nrow = 2, ncol = 2)) +
theme_classic(base_size = 15) +
theme(legend.position = "bottom",
legend.background = element_rect(color = "black")) +
guides(color = FALSE)
Kaplan_TrT_CYBA
survdiff(Surv(Weeks.Alive, Status) ~ Treatment, data = Fern2.CYBA)
## Call:
## survdiff(formula = Surv(Weeks.Alive, Status) ~ Treatment, data = Fern2.CYBA)
##
## n=77, 158 observations deleted due to missingness.
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## Treatment=No Competition 13 10 7.96 0.523 0.6438
## Treatment=Only Belowground Competition 12 9 8.57 0.022 0.0275
## Treatment=Only Aboveground Competition 22 14 19.39 1.497 2.3932
## Treatment=Full Competition 30 25 22.09 0.384 0.6570
##
## Chisq= 2.6 on 3 degrees of freedom, p= 0.5
png("~/Downloads/KaplanTreatment_CYBA.png", width = 165, height = 120, units = 'mm', res = 300)
Kaplan_TrT_CYBA
dev.off()
## quartz_off_screen
## 2
fi_trt_fit_CYBA_AG <- survfit(Surv(Weeks.Alive, Status) ~ Aboveground, data=Fern2.CYBA)
Kaplan_TrT_CYBA_AG <- autoplot(fi_trt_fit_CYBA_AG) +
ylab("Probability of Survival") +
xlab("Weeks Since Germination") +
geom_label(aes(x = 1, y = 0,
label = paste0("p = 0.5")),
vjust = -9.6, hjust = -5.825,
size = 5, fill = "white", color = "black") +
scale_color_manual(values = c("black", "black", "black", "black"),
guide = guide_legend(direction = "horizontal",
title = NULL,
nrow = 2, ncol = 2)) +
scale_fill_manual(values = c("hotpink", "#1F78B4", "limegreen", "#FFD700"),
guide = guide_legend(direction = "horizontal",
title = NULL,
nrow = 2, ncol = 2)) +
theme_classic(base_size = 15) +
theme(legend.position = "bottom",
legend.background = element_rect(color = "black")) +
guides(color = FALSE)
Kaplan_TrT_CYBA_AG
survdiff(Surv(Weeks.Alive, Status) ~ Aboveground, data = Fern2.CYBA)
## Call:
## survdiff(formula = Surv(Weeks.Alive, Status) ~ Aboveground, data = Fern2.CYBA)
##
## n=77, 158 observations deleted due to missingness.
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## Aboveground=Covered 52 39 41.5 0.148 0.551
## Aboveground=Exposed 25 19 16.5 0.371 0.551
##
## Chisq= 0.6 on 1 degrees of freedom, p= 0.5
png("~/Downloads/KaplanAG_CYBA.png", width = 165, height = 120, units = 'mm', res = 300)
Kaplan_TrT_CYBA_AG
dev.off()
## quartz_off_screen
## 2
fi_trt_fit_CYBA_BG <- survfit(Surv(Weeks.Alive, Status) ~ Belowground, data=Fern2.CYBA)
Kaplan_TrT_CYBA_BG <- autoplot(fi_trt_fit_CYBA_BG) +
ylab("Probability of Survival") +
xlab("Weeks Since Germination") +
geom_label(aes(x = 1, y = 0,
label = paste0("p = 0.4")),
vjust = -9.6, hjust = -5.825,
size = 5, fill = "white", color = "black") +
scale_color_manual(values = c("black", "black", "black", "black"),
guide = guide_legend(direction = "horizontal",
title = NULL,
nrow = 2, ncol = 2)) +
scale_fill_manual(values = c("hotpink", "#1F78B4", "limegreen", "#FFD700"),
guide = guide_legend(direction = "horizontal",
title = NULL,
nrow = 2, ncol = 2)) +
theme_classic(base_size = 15) +
theme(legend.position = "bottom",
legend.background = element_rect(color = "black")) +
guides(color = FALSE)
Kaplan_TrT_CYBA_BG
survdiff(Surv(Weeks.Alive, Status) ~ Belowground, data = Fern2.CYBA)
## Call:
## survdiff(formula = Surv(Weeks.Alive, Status) ~ Belowground, data = Fern2.CYBA)
##
## n=77, 158 observations deleted due to missingness.
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## Belowground=Protected 35 24 27.3 0.410 0.824
## Belowground=Unprotected 42 34 30.7 0.365 0.824
##
## Chisq= 0.8 on 1 degrees of freedom, p= 0.4
png("~/Downloads/KaplanBG_CYBA.png", width = 165, height = 120, units = 'mm', res = 300)
Kaplan_TrT_CYBA_BG
dev.off()
## quartz_off_screen
## 2
fi_trt_fit <- survfit(Surv(Weeks.Alive, Status) ~ Treatment, data=Fern2)
Kaplan_TrT <- autoplot(fi_trt_fit) +
ylab("Probability of Survival") +
xlab("Weeks Since Germination") +
geom_label(aes(x = 1, y = 0,
label = paste0("p = 0.3")),
vjust = -9.6, hjust = -5.825,
size = 5, fill = "white", color = "black") +
scale_color_manual(values = c("black", "black", "black", "black"),
guide = guide_legend(direction = "horizontal",
title = NULL,
nrow = 2, ncol = 2)) +
scale_fill_manual(values = c("hotpink", "#1F78B4", "limegreen", "#FFD700"),
guide = guide_legend(direction = "horizontal",
title = NULL,
nrow = 2, ncol = 2)) +
theme_classic(base_size = 15) +
theme(legend.position = "bottom",
legend.background = element_rect(color = "black")) +
guides(color = FALSE)
Kaplan_TrT
survdiff(Surv(Weeks.Alive, Status) ~ Treatment, data = Fern2)
## Call:
## survdiff(formula = Surv(Weeks.Alive, Status) ~ Treatment, data = Fern2)
##
## n=122, 348 observations deleted due to missingness.
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## Treatment=No Competition 26 16 14.9 0.0754 0.0986
## Treatment=Only Belowground Competition 22 12 13.7 0.2216 0.2837
## Treatment=Only Aboveground Competition 32 16 21.5 1.3863 2.0435
## Treatment=Full Competition 42 29 22.9 1.6479 2.5037
##
## Chisq= 3.5 on 3 degrees of freedom, p= 0.3
png("~/Downloads/KaplanTreatment.png", width = 165, height = 120, units = 'mm', res = 300)
Kaplan_TrT
dev.off()
## quartz_off_screen
## 2
fi_spp_fit <- survfit(Surv(Weeks.Alive, Status) ~ Species, data=Fern2)
Kaplan_Spp <- autoplot(fi_spp_fit) +
ylab("Probability of Survival") +
xlab("Weeks Since Germination") +
geom_label(aes(x = 1, y = 0,
label = paste0("p = 4e-05")),
vjust = -10.75, hjust = -4.175,
size = 5, fill = "white", color = "black") +
scale_color_manual(values = c("black", "black", "black", "black"),
guide = guide_legend(direction = "horizontal",
title = NULL,
nrow = 2, ncol = 2)) +
scale_fill_manual(values = c("#bef7f1", "#93A02C"),
guide = guide_legend(direction = "horizontal",
title = NULL,
nrow = 1, ncol = 2)) +
theme_classic(base_size = 15) +
theme(legend.position = "bottom",
legend.background = element_rect(color = "black"),
legend.text = element_text(face = "italic")) + # This line italicizes the legend labels) +
guides(color = FALSE)
Kaplan_Spp
survdiff(Surv(Weeks.Alive, Status) ~ Species, data = Fern2)
## Call:
## survdiff(formula = Surv(Weeks.Alive, Status) ~ Species, data = Fern2)
##
## n=122, 348 observations deleted due to missingness.
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## Species=Cymbopetalum baillonii 77 58 41.2 6.89 16.7
## Species=Sapindus saponaria 45 15 31.8 8.90 16.7
##
## Chisq= 16.7 on 1 degrees of freedom, p= 4e-05
png("~/Downloads/KaplanSpecies.png", width = 165, height = 120, units = 'mm', res = 300)
Kaplan_Spp
dev.off()
## quartz_off_screen
## 2
Seedling Establishment between
Species?glm.nb1_spp <- glm.nb(data = T7B, Count ~ Species + Aboveground + Belowground)
glm_TMB_spp <- glmmTMB(Count ~ Species + Aboveground + Belowground + (1 | Plot),
family = nbinom2(),
data = T7B)
Anova(glm.nb1_spp)
## Analysis of Deviance Table (Type II tests)
##
## Response: Count
## LR Chisq Df Pr(>Chisq)
## Species 8.231 1 0.00412 **
## Aboveground 2.733 1 0.09827 .
## Belowground 0.020 1 0.88804
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(glm_TMB_spp)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Count
## Chisq Df Pr(>Chisq)
## Species 7.954 1 0.0048 **
## Aboveground 2.678 1 0.1018
## Belowground 0.020 1 0.8884
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(glm_TMB_spp)
## Family: nbinom2 ( log )
## Formula: Count ~ Species + Aboveground + Belowground + (1 | Plot)
## Data: T7B
##
## AIC BIC logLik deviance df.resid
## 192.9 204.1 -90.4 180.9 42
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 1.29e-09 3.59e-05
## Number of obs: 48, groups: Plot, 6
##
## Dispersion parameter for nbinom2 family (): 35.5
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.3540 0.1684 8.04 8.9e-16 ***
## SpeciesSapindus saponaria -0.5317 0.1885 -2.82 0.0048 **
## AbovegroundExposed -0.3034 0.1854 -1.64 0.1018
## BelowgroundUnprotected 0.0257 0.1831 0.14 0.8884
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Answer: Yes, there is a significant difference
between Species on Seedling Establishment
(X2 = 8.231, p = 0.00412). Version with RE:
(X2 = 7.954, p = 0.0048).
# Given values
var_random_effect <- 1.29e-09 # Variance of the random effect for Plot
residual_variance <- 3.29 # Approximated residual variance
# Calculate ICC
ICC <- var_random_effect / (var_random_effect + residual_variance)
# Convert to percentage
ICC_percent <- ICC * 100
ICC_percent
## [1] 3.92097e-08
Belowground Treatment and
Aboveground Treatment on
Seedling Establishment?T7B.SASA <- droplevels(subset(T7B, T7B$Species == "Sapindus saponaria"))
glm.nb1_SASA <- glm.nb(data = T7B.SASA, Count ~ Aboveground * Belowground)
glm_TMB_SASA <- glmmTMB(Count ~ Aboveground * Belowground + (1 | Plot),
family = nbinom2(),
data = T7B.SASA)
Anova(glm.nb1_SASA)
## Analysis of Deviance Table (Type II tests)
##
## Response: Count
## LR Chisq Df Pr(>Chisq)
## Aboveground 0.3337 1 0.563
## Belowground 0.0834 1 0.773
## Aboveground:Belowground 0.7165 1 0.397
Anova(glm_TMB_SASA)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Count
## Chisq Df Pr(>Chisq)
## Aboveground 0.318 1 0.573
## Belowground 0.080 1 0.777
## Aboveground:Belowground 0.711 1 0.399
Answer: No, there was no significant interaction
between Belowground Treatment and
Aboveground Treatment on
Seedling Establishment for Sapindus saponaria
(X2 = 0.7165, p = 0.397). Version with RE: (X2 =
0.711, p = 0.399).
Seedling Establishment between
Belowground Treatments for Sapindus
saponaria?Seedling Establishment between
Aboveground Treatments for Sapindus
saponaria?glm.nb1_SASA2 <- glm.nb(data = T7B.SASA, Count ~ Aboveground + Belowground)
glm.TMB_SASA2 <- glmmTMB(Count ~ Aboveground + Belowground + (1 | Plot),
family = nbinom2(),
data = T7B.SASA)
Anova(glm.nb1_SASA2)
## Analysis of Deviance Table (Type II tests)
##
## Response: Count
## LR Chisq Df Pr(>Chisq)
## Aboveground 0.3337 1 0.563
## Belowground 0.0834 1 0.773
Anova(glm.TMB_SASA2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Count
## Chisq Df Pr(>Chisq)
## Aboveground 0.333 1 0.564
## Belowground 0.083 1 0.773
summary(glm.TMB_SASA2)
## Family: nbinom2 ( log )
## Formula: Count ~ Aboveground + Belowground + (1 | Plot)
## Data: T7B.SASA
##
## AIC BIC logLik deviance df.resid
## 86.2 92.1 -38.1 76.2 19
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 7.93e-10 2.82e-05
## Number of obs: 24, groups: Plot, 6
##
## Dispersion parameter for nbinom2 family (): 8.83e+07
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.6470 0.2542 2.54 0.011 *
## AbovegroundExposed 0.1671 0.2897 0.58 0.564
## BelowgroundUnprotected -0.0834 0.2889 -0.29 0.773
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Answer:
No, there is not a significant difference between
Aboveground Treatments on
Seedling Establishment for Sapindus saponaria
(X2 = 0.3337, p = 0.563). Version with RE:
(X2 = 0.333, p = 0.564).
No, there is a no significant difference between
Belowground Treatments on
Seedling Establishment for Sapindus saponaria
(X2 = 0.0834, p = 0.773). Version with RE:
(X2 = 0.083, p = 0.773).
Belowground Treatment and
Aboveground Treatment on
Seedling Establishment?T7B.CYBA <- droplevels(subset(T7B, T7B$Species == "Cymbopetalum baillonii"))
glm.nb1_CYBA <- glm.nb(data = T7B.CYBA, Count ~ Aboveground * Belowground)
glm.TMB_CYBA <- glmmTMB(Count ~ Aboveground * Belowground + (1 | Plot),
family = nbinom2(),
data = T7B.CYBA)
Anova(glm.nb1_CYBA)
## Analysis of Deviance Table (Type II tests)
##
## Response: Count
## LR Chisq Df Pr(>Chisq)
## Aboveground 6.795 1 0.00914 **
## Belowground 0.170 1 0.67988
## Aboveground:Belowground 2.116 1 0.14577
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(glm.TMB_CYBA)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Count
## Chisq Df Pr(>Chisq)
## Aboveground 6.024 1 0.0141 *
## Belowground 0.170 1 0.6801
## Aboveground:Belowground 2.088 1 0.1485
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(glm.TMB_CYBA)
## Family: nbinom2 ( log )
## Formula: Count ~ Aboveground * Belowground + (1 | Plot)
## Data: T7B.CYBA
##
## AIC BIC logLik deviance df.resid
## 110.0 117.1 -49.0 98.0 18
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 0.000353 0.0188
## Number of obs: 24, groups: Plot, 6
##
## Dispersion parameter for nbinom2 family (): 74.2
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.299 0.220 5.91 3.4e-09
## AbovegroundExposed -0.258 0.330 -0.78 0.43
## BelowgroundUnprotected 0.343 0.287 1.20 0.23
## AbovegroundExposed:BelowgroundUnprotected -0.691 0.478 -1.44 0.15
##
## (Intercept) ***
## AbovegroundExposed
## BelowgroundUnprotected
## AbovegroundExposed:BelowgroundUnprotected
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Answer: No, there was no significant interaction
between Belowground Treatment and
Aboveground Treatment on
Seedling Establishment for Cymbopetalum baillonii
(X2 = 2.116, p = 0.14577). Version with RE:
(X2 = 2.088, p = 0.1485)
Seedling Establishment between
Belowground Treatments for Cymbopetalum
baillonii?Seedling Establishment between
Aboveground Treatments for Cymbopetalum
baillonii?glm.nb1_CYBA2 <- glm.nb(data = T7B.CYBA, Count ~ Aboveground + Belowground)
glm.TMB_CYBA2 <- glmmTMB(Count ~ Aboveground + Belowground + (1 | Plot),
family = nbinom2(),
data = T7B.CYBA)
Anova(glm.nb1_CYBA2)
## Analysis of Deviance Table (Type II tests)
##
## Response: Count
## LR Chisq Df Pr(>Chisq)
## Aboveground 6.111 1 0.0134 *
## Belowground 0.125 1 0.7231
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(glm.TMB_CYBA2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Count
## Chisq Df Pr(>Chisq)
## Aboveground 5.899 1 0.0151 *
## Belowground 0.124 1 0.7249
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Answer:
Yes, there is not a significant difference between
Aboveground Treatments on
Seedling Establishment for Cymbopetalum baillonii
(X2 = 6.111, p = 0.0134). Version with RE:
(X2 = 5.899, p = 0.0151)
No, there is a no significant difference between
Belowground Treatments on
Seedling Establishment for Cymbopetalum baillonii
(X2 = 0.125, p = 0.7231). Version with RE:
(X2 = 0.124, p = 0.7249)
Maximum Leaf Count between
Species?zinb_3 <- zeroinfl(MaxLeaves ~ Species + Aboveground + Belowground | 1, data = FernExp.G2, dist = "negbin")
Anova(zinb_3)
## Analysis of Deviance Table (Type II tests)
##
## Response: MaxLeaves
## Df Chisq Pr(>Chisq)
## Species 1 4.298 0.0382 *
## Aboveground 1 3.896 0.0484 *
## Belowground 1 1.009 0.3151
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
zinb_3_glmmTMB <- glmmTMB(MaxLeaves ~ Species + Aboveground + Belowground + (1 | Plot),
ziformula = ~ 1,
family = nbinom2(),
data = FernExp.G2)
Anova(zinb_3_glmmTMB)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: MaxLeaves
## Chisq Df Pr(>Chisq)
## Species 4.354 1 0.0369 *
## Aboveground 4.356 1 0.0369 *
## Belowground 1.325 1 0.2497
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(zinb_3_glmmTMB)
## Family: nbinom2 ( log )
## Formula: MaxLeaves ~ Species + Aboveground + Belowground + (1 | Plot)
## Zero inflation: ~1
## Data: FernExp.G2
##
## AIC BIC logLik deviance df.resid
## 704.2 724.0 -345.1 690.2 118
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 0.0188 0.137
## Number of obs: 125, groups: Plot, 6
##
## Dispersion parameter for nbinom2 family (): 1.72
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.643 0.175 9.38 <2e-16 ***
## SpeciesSapindus saponaria 0.389 0.186 2.09 0.037 *
## AbovegroundExposed 0.413 0.198 2.09 0.037 *
## BelowgroundUnprotected -0.210 0.182 -1.15 0.250
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.735 0.366 -4.74 2.2e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Answer: Yes, there is a significant difference
between Species on Leaf Count (X2 =
4.298, p = 0.0382). Version with RE: (X2 =
4.354, p = 0.0369)
To calculate the Intraclass Correlation Coefficient (ICC) for this model, we need to know both the variance of the random effect and the residual variance. In a generalized linear mixed model (GLMM) with a negative binomial distribution, the residual variance isn’t directly provided in the model summary. However, we can estimate it using the dispersion parameter (theta) and the link function (log in this case).
Residual Variance= π2/3 ≈ 3.29
This approximation works for the logit link in binomial models, but is commonly adapted for GLMMs with log links in the absence of a direct residual variance estimate.
# Given values
var_random_effect <- 0.0188 # Variance of the random effect for Plot
residual_variance <- 3.29 # Approximated residual variance
# Calculate ICC
ICC <- var_random_effect / (var_random_effect + residual_variance)
ICC*100
## [1] 0.568182
Answer: The ICC is approximately 0.0057, or 0.57%.
This suggests that only about 0.57% of the total variation in
MaxLeaves is attributable to differences between plots.
This further supports the notion that plot-to-plot variability is
relatively small compared to the overall variability in the data.
Belowground Treatment and
Aboveground Treatment on
Maximum Leaf Count?Anova(zinb_SASA)
## Analysis of Deviance Table (Type II tests)
##
## Response: MaxLeaves
## Df Chisq Pr(>Chisq)
## Aboveground 1 10.327 0.00131 **
## Belowground 1 0.421 0.51661
## Aboveground:Belowground 1 0.463 0.49617
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(zinb_SASA_GLMMTMB)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: MaxLeaves
## Chisq Df Pr(>Chisq)
## Aboveground 11.141 1 0.000844 ***
## Belowground 0.236 1 0.627316
## Aboveground:Belowground 0.342 1 0.558950
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(zinb_SASA_GLMMTMB)
## Family: nbinom2 ( log )
## Formula: MaxLeaves ~ Aboveground * Belowground + (1 | Plot)
## Zero inflation: ~1
## Data: FernExp.G.SASA
##
## AIC BIC logLik deviance df.resid
## 293.1 306.0 -139.5 279.1 40
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 0.0379 0.195
## Number of obs: 47, groups: Plot, 6
##
## Dispersion parameter for nbinom2 family (): 4.07
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.693 0.230 7.37 1.8e-13
## AbovegroundExposed 0.771 0.274 2.81 0.0049
## BelowgroundUnprotected 0.238 0.314 0.76 0.4474
## AbovegroundExposed:BelowgroundUnprotected -0.233 0.398 -0.58 0.5589
##
## (Intercept) ***
## AbovegroundExposed **
## BelowgroundUnprotected
## AbovegroundExposed:BelowgroundUnprotected
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.004 0.473 -4.23 2.3e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Answer: No, there is no significant interaction
between Belowground Treatment and
Aboveground Treatment on Leaf Count for
Sapindus saponaria (X2 = 0.463, p = 0.496). Version
with RE: (X2 = 0.342, p = 0.559)
Leaf Count between
Belowground Treatments for Sapindus
saponaria?Leaf Count between
Aboveground Treatments for Sapindus
saponaria?zinb_SASA2 <- zeroinfl(MaxLeaves ~ Aboveground + Belowground | 1, data = FernExp.G.SASA, dist = "negbin")
zinb_SASA2_GLMMTMB <- glmmTMB(MaxLeaves ~ Aboveground + Belowground + (1 | Plot),
ziformula = ~ 1,
family = nbinom2(),
data = FernExp.G.SASA)
Anova(zinb_SASA2)
## Analysis of Deviance Table (Type II tests)
##
## Response: MaxLeaves
## Df Chisq Pr(>Chisq)
## Aboveground 1 10.101 0.00148 **
## Belowground 1 0.411 0.52125
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(zinb_SASA2_GLMMTMB)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: MaxLeaves
## Chisq Df Pr(>Chisq)
## Aboveground 11.028 1 0.000897 ***
## Belowground 0.233 1 0.629195
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Answer:
Yes, there is a significant difference between
Aboveground Treatments on Leaf Count for
Sapindus saponaria (X2 = 10.101, p =
0.0015). Version with RE: (X2 = 11.028, p =
0.0009)
No, there is a no significant difference between
Belowground Treatments on Leaf Count for
Sapindus saponaria (X2 = 0.411, p =
0.521). Version with RE: (X2 = 0.233, p =
0.629)
Belowground Treatment and
Aboveground Treatment on
Maximum Leaf Count?Anova(zinb_CYBA)
## Analysis of Deviance Table (Type II tests)
##
## Response: MaxLeaves
## Df Chisq Pr(>Chisq)
## Aboveground 1 0.189 0.664
## Belowground 1 1.778 0.182
## Aboveground:Belowground 1 0.026 0.873
Anova(zinb_CYBA_GLMMTMB)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: MaxLeaves
## Chisq Df Pr(>Chisq)
## Aboveground 0.189 1 0.664
## Belowground 1.778 1 0.182
## Aboveground:Belowground 0.026 1 0.873
Answer: No, there is no significant interaction
between Belowground Treatment and
Aboveground Treatment on Leaf Count for
Cymbopetalum baillonii (X2 = 0.026, p =
0.873). Version with RE: (X2 = 0.026, p =
0.873).
Leaf Count between
Belowground Treatments for Cymbopetalum
baillonii?Leaf Count between
Aboveground Treatments for Cymbopetalum
baillonii?zinb_CYBA2 <- zeroinfl(MaxLeaves ~ Aboveground + Belowground | 1, data = FernExp.G.CYBA, dist = "negbin")
zinb_CYBA2_GLMMTMB <- glmmTMB(MaxLeaves ~ Aboveground + Belowground + (1 | Plot),
ziformula = ~ 1,
family = nbinom2(),
data = FernExp.G.CYBA)
Anova(zinb_CYBA2)
## Analysis of Deviance Table (Type II tests)
##
## Response: MaxLeaves
## Df Chisq Pr(>Chisq)
## Aboveground 1 0.191 0.662
## Belowground 1 1.794 0.180
Anova(zinb_CYBA2_GLMMTMB)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: MaxLeaves
## Chisq Df Pr(>Chisq)
## Aboveground 0.191 1 0.662
## Belowground 1.794 1 0.180
summary(zinb_CYBA2_GLMMTMB)
## Family: nbinom2 ( log )
## Formula: MaxLeaves ~ Aboveground + Belowground + (1 | Plot)
## Zero inflation: ~1
## Data: FernExp.G.CYBA
##
## AIC BIC logLik deviance df.resid
## 410.3 424.4 -199.1 398.3 72
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 2.54e-09 5.04e-05
## Number of obs: 78, groups: Plot, 6
##
## Dispersion parameter for nbinom2 family (): 1.12
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.799 0.242 7.44 1e-13 ***
## AbovegroundExposed 0.125 0.287 0.44 0.66
## BelowgroundUnprotected -0.364 0.272 -1.34 0.18
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.742 0.683 -2.55 0.011 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Answer:
No, there is a significant difference between
Aboveground Treatments on Leaf Count for
Cymbopetalum baillonii (X2 = 0.191, p =
0.662). Version with RE: (X2 = 0.191, p =
0.662).
No, there is a no significant difference between
Belowground Treatments on Leaf Count for
Cymbopetalum baillonii (X2 = 1.794, p =
0.180). Version with RE: (X2 = 1.794, p =
0.180).
WetMassExp.SASA <- droplevels(subset(WetMassExp, WetMassExp$Species == "Sapindus saponaria"))
WetMassExp.CYBA <- droplevels(subset(WetMassExp, WetMassExp$Species == "Cymbopetalum baillonii"))
glmHeight.SASA <- glm(data = WetMassExp.SASA, Height..cm. ~ Aboveground * Belowground,
family = Gamma(link = "log"))
glmHeight.CYBA <- glm(data = WetMassExp.CYBA, Height..cm. ~ Aboveground * Belowground,
family = Gamma(link = "log"))
glmHeight.SASA_GLMMTMB <- glmmTMB(Height..cm. ~ Aboveground * Belowground + (1 | Plot),
family = Gamma(link = "log"),
data = WetMassExp.SASA)
glmHeight.CYBA_GLMMTMB <- glmmTMB(Height..cm. ~ Aboveground * Belowground + (1 | Plot),
family = Gamma(link = "log"),
data = WetMassExp.CYBA)
Seedling Height between
species?glmHeight_spp <- glm(data = WetMassExp, Height..cm. ~ Species + Aboveground + Belowground,
family = Gamma(link = "log"))
glmHeight_spp_GLMMTMB <- glmmTMB(Height..cm. ~ Species + Aboveground + Belowground + (1 | Plot),
family = Gamma(link = "log"),
data = WetMassExp)
Anova(glmHeight_spp)
## Analysis of Deviance Table (Type II tests)
##
## Response: Height..cm.
## LR Chisq Df Pr(>Chisq)
## Species 0.008 1 0.93091
## Aboveground 7.852 1 0.00508 **
## Belowground 0.487 1 0.48542
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(glmHeight_spp_GLMMTMB)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Height..cm.
## Chisq Df Pr(>Chisq)
## Species 0.006 1 0.937
## Aboveground 18.529 1 1.67e-05 ***
## Belowground 0.077 1 0.781
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(glmHeight_spp_GLMMTMB)
## Family: Gamma ( log )
## Formula:
## Height..cm. ~ Species + Aboveground + Belowground + (1 | Plot)
## Data: WetMassExp
##
## AIC BIC logLik deviance df.resid
## 378.8 390.3 -183.4 366.8 44
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 0.0948 0.308
## Number of obs: 50, groups: Plot, 6
##
## Dispersion estimate for Gamma family (sigma^2): 0.216
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.7243 0.1849 14.73 < 2e-16 ***
## SpeciesSapindus saponaria 0.0115 0.1456 0.08 0.94
## AbovegroundExposed 0.6547 0.1521 4.30 1.7e-05 ***
## BelowgroundUnprotected -0.0407 0.1464 -0.28 0.78
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Answer: No, there is not a significant difference
between Species on Seedling Height
(X2 = 0.487, p = 0.93091). Version with RE:
(X2 = 0.006, p = 0.937)
Belowground Treatment and
Aboveground Treatment on Seedling Height?Anova(glmHeight.SASA)
## Analysis of Deviance Table (Type II tests)
##
## Response: Height..cm.
## LR Chisq Df Pr(>Chisq)
## Aboveground 9.100 1 0.00256 **
## Belowground 0.492 1 0.48300
## Aboveground:Belowground 0.006 1 0.93764
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(glmHeight.SASA_GLMMTMB)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Height..cm.
## Chisq Df Pr(>Chisq)
## Aboveground 18.441 1 1.75e-05 ***
## Belowground 0.024 1 0.876
## Aboveground:Belowground 0.069 1 0.792
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Answer: No, there was no significant interaction
between Belowground Treatment and
Aboveground Treatment on Seedling Height for
Sapindus saponaria (X2 = 0.006, p =
0.938). Version with RE: (X2 = 0.069, p =
0.792).
Seedling Height
between Belowground Treatments for Sapindus
saponaria?Seedling Height
between Aboveground Treatments for Sapindus
saponaria?glmHeight.SASA2 <- glm(data = WetMassExp.SASA, Height..cm. ~ Aboveground + Belowground,
family = Gamma(link = "log"))
glmHeight.SASA_GLMMTMB2 <- glmmTMB(Height..cm. ~ Aboveground + Belowground + (1 | Plot),
family = Gamma(link = "log"),
data = WetMassExp.SASA)
Anova(glmHeight.SASA2)
## Analysis of Deviance Table (Type II tests)
##
## Response: Height..cm.
## LR Chisq Df Pr(>Chisq)
## Aboveground 9.527 1 0.00202 **
## Belowground 0.515 1 0.47291
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(glmHeight.SASA_GLMMTMB2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Height..cm.
## Chisq Df Pr(>Chisq)
## Aboveground 18.357 1 1.83e-05 ***
## Belowground 0.024 1 0.876
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(glmHeight.SASA_GLMMTMB2)
## Family: Gamma ( log )
## Formula: Height..cm. ~ Aboveground + Belowground + (1 | Plot)
## Data: WetMassExp.SASA
##
## AIC BIC logLik deviance df.resid
## 241.5 248.7 -115.7 231.5 26
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 0.129 0.359
## Number of obs: 31, groups: Plot, 6
##
## Dispersion estimate for Gamma family (sigma^2): 0.227
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.6533 0.2143 12.38 < 2e-16 ***
## AbovegroundExposed 0.7701 0.1797 4.28 1.8e-05 ***
## BelowgroundUnprotected 0.0294 0.1877 0.16 0.88
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Answer:
Yes, there is a significant difference between
Aboveground Treatments on Seedling Height for
Sapindus saponaria (X2 = 9.527, p =
0.00202). Version with RE: (X2 = 18.357, p =
1.83e-05).
No, there is a no significant difference between
Belowground Treatments on Seedling Height for
Sapindus saponaria (X2 = 0.515, p =
0.47291). Version with RE: (X2 = 0.024, p =
0.876).
Belowground Treatment and
Aboveground Treatment on on
Seedling Height?Anova(glmHeight.CYBA)
## Analysis of Deviance Table (Type II tests)
##
## Response: Height..cm.
## LR Chisq Df Pr(>Chisq)
## Aboveground 0.13673 1 0.712
## Belowground 0.19647 1 0.658
## Aboveground:Belowground 0.01569 1 0.900
Anova(glmHeight.CYBA_GLMMTMB)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Height..cm.
## Chisq Df Pr(>Chisq)
## Aboveground 0.147 1 0.701
## Belowground 0.214 1 0.644
## Aboveground:Belowground 0.017 1 0.896
Answer: No, there is no significant interaction
between Belowground Treatment and
Aboveground Treatment on Seedling Height
(X2 = 0.0157, p = 0.900). Version with RE:
(X2 = 0.017, p = 0.896).
Seedling Height
between Belowground Treatments for Cymbopetalum
baillonii?Seedling Height
between Aboveground Treatments for Cymbopetalum
baillonii?glmHeight.CYBA2 <- glm(data = WetMassExp.CYBA, Height..cm. ~ Aboveground + Belowground,
family = Gamma(link = "log"))
glmHeight.CYBA_GLMMTMB2 <- glmmTMB(Height..cm. ~ Aboveground + Belowground + (1 | Plot),
family = Gamma(link = "log"),
data = WetMassExp.CYBA)
Anova(glmHeight.CYBA2)
## Analysis of Deviance Table (Type II tests)
##
## Response: Height..cm.
## LR Chisq Df Pr(>Chisq)
## Aboveground 0.1453 1 0.703
## Belowground 0.2088 1 0.648
Anova(glmHeight.CYBA_GLMMTMB2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Height..cm.
## Chisq Df Pr(>Chisq)
## Aboveground 0.147 1 0.701
## Belowground 0.213 1 0.644
summary(glmHeight.CYBA_GLMMTMB2)
## Family: Gamma ( log )
## Formula: Height..cm. ~ Aboveground + Belowground + (1 | Plot)
## Data: WetMassExp.CYBA
##
## AIC BIC logLik deviance df.resid
## 138.0 142.7 -64.0 128.0 14
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 1.6e-10 1.27e-05
## Number of obs: 19, groups: Plot, 5
##
## Dispersion estimate for Gamma family (sigma^2): 0.144
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.9068 0.1320 22.01 <2e-16 ***
## AbovegroundExposed 0.0719 0.1874 0.38 0.70
## BelowgroundUnprotected 0.0805 0.1744 0.46 0.64
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Answer:
No, there is a significant difference between
Aboveground Treatments on Leaf Count for
Cymbopetalum baillonii (X2 = 0.1453, p =
0.703). Version with RE: (X2 = 0.147, p =
0.701).
No, there is a no significant difference between
Belowground Treatments on Leaf Count for
Cymbopetalum baillonii (X2 = 0.2088, p =
0.648). Version with RE: (X2 = 0.213, p =
0.644).
glmG.Diam.SASA <- glm(data = WetMassExp.SASA, Diameter.Base..cm. ~ Aboveground * Belowground,
family = Gamma(link = "log"))
glmG.Diam.CYBA <- glm(data = WetMassExp.CYBA, Diameter.Base..cm. ~ Aboveground * Belowground,
family = Gamma(link = "log"))
glmG.Diam.SASA_glmmTMB<- glmmTMB(Diameter.Base..cm. ~ Aboveground * Belowground + (1 | Plot),
family = Gamma(link = "log"),
data = WetMassExp.SASA)
glmG.Diam.CYBA_glmmTMB<- glmmTMB(Diameter.Base..cm. ~ Aboveground * Belowground + (1 | Plot),
family = Gamma(link = "log"),
data = WetMassExp.CYBA)
WetMassExp.C.SASA <- droplevels(subset(WetMassExp.C, WetMassExp.C$Species == "Sapindus saponaria"))
glmG.Diam.C.SASA <- glm(data = WetMassExp.C.SASA, Diameter.Base..cm. ~ Aboveground * Belowground,
family = Gamma(link = "log"))
glmG.Diam.C.SASA_glmmTMB <- glmmTMB(Diameter.Base..cm. ~ Aboveground * Belowground + (1 | Plot),
family = Gamma(link = "log"),
data = WetMassExp.C.SASA)
Belowground Treatment and
Aboveground Treatment on Base Diameter for
Sapindus saponaria? What about when we exclude the outlier? In
the Only Belowground Competition treatment, we had one very
large Sapindus.Base Diameter between
species?glmG.Diam2 <- glm(data = WetMassExp, Diameter.Base..cm. ~ Species + Aboveground + Belowground,
family = Gamma(link = "log"))
glmG.Diam2_glmmTMB<- glmmTMB(Diameter.Base..cm. ~ Species + Aboveground + Belowground + (1 | Plot),
family = Gamma(link = "log"),
data = WetMassExp)
Anova(glmG.Diam2)
## Analysis of Deviance Table (Type II tests)
##
## Response: Diameter.Base..cm.
## LR Chisq Df Pr(>Chisq)
## Species 0.361 1 0.547703
## Aboveground 11.429 1 0.000723 ***
## Belowground 2.561 1 0.109520
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(glmG.Diam2_glmmTMB)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Diameter.Base..cm.
## Chisq Df Pr(>Chisq)
## Species 1.380 1 0.2400
## Aboveground 35.878 1 2.1e-09 ***
## Belowground 4.071 1 0.0436 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(glmG.Diam2_glmmTMB)
## Family: Gamma ( log )
## Formula:
## Diameter.Base..cm. ~ Species + Aboveground + Belowground + (1 | Plot)
## Data: WetMassExp
##
## AIC BIC logLik deviance df.resid
## 23.5 34.9 -5.7 11.5 44
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 0.132 0.363
## Number of obs: 50, groups: Plot, 6
##
## Dispersion estimate for Gamma family (sigma^2): 0.277
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.388 0.212 -6.54 6.3e-11 ***
## SpeciesSapindus saponaria 0.200 0.170 1.17 0.240
## AbovegroundExposed 1.073 0.179 5.99 2.1e-09 ***
## BelowgroundUnprotected 0.322 0.160 2.02 0.044 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Answer: No, there is not a significant difference
between Species on Base Diameter
(X2 = 0.361, p = 0.548). Version with RE:
(X2 = 1.38, p = 0.24)
Anova(glmG.Diam.SASA)
## Analysis of Deviance Table (Type II tests)
##
## Response: Diameter.Base..cm.
## LR Chisq Df Pr(>Chisq)
## Aboveground 19.664 1 9.23e-06 ***
## Belowground 4.898 1 0.0269 *
## Aboveground:Belowground 2.917 1 0.0876 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(glmG.Diam.C.SASA)
## Analysis of Deviance Table (Type II tests)
##
## Response: Diameter.Base..cm.
## LR Chisq Df Pr(>Chisq)
## Aboveground 21.586 1 3.38e-06 ***
## Belowground 0.947 1 0.330
## Aboveground:Belowground 0.113 1 0.737
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(glmG.Diam.SASA_glmmTMB)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Diameter.Base..cm.
## Chisq Df Pr(>Chisq)
## Aboveground 32.460 1 1.22e-08 ***
## Belowground 4.579 1 0.0324 *
## Aboveground:Belowground 1.458 1 0.2272
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(glmG.Diam.C.SASA_glmmTMB)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Diameter.Base..cm.
## Chisq Df Pr(>Chisq)
## Aboveground 43.872 1 3.51e-11 ***
## Belowground 0.478 1 0.489
## Aboveground:Belowground 0.098 1 0.754
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(glmG.Diam.C.SASA_glmmTMB)
## Family: Gamma ( log )
## Formula: Diameter.Base..cm. ~ Aboveground * Belowground + (1 | Plot)
## Data: WetMassExp.C.SASA
##
## AIC BIC logLik deviance df.resid
## -0.8 7.6 6.4 -12.8 24
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 0.0498 0.223
## Number of obs: 30, groups: Plot, 6
##
## Dispersion estimate for Gamma family (sigma^2): 0.151
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.2188 0.1704 -7.15 8.5e-13
## AbovegroundExposed 0.9382 0.1985 4.73 2.3e-06
## BelowgroundUnprotected 0.0624 0.2120 0.29 0.77
## AbovegroundExposed:BelowgroundUnprotected 0.0915 0.2920 0.31 0.75
##
## (Intercept) ***
## AbovegroundExposed ***
## BelowgroundUnprotected
## AbovegroundExposed:BelowgroundUnprotected
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Answer: No, there is no significant interaction
between Belowground Treatment and
Aboveground Treatment on Base Diameter for
Sapindus saponaria (X2 = 2.917, p =
0.0876). When we exclude this outlier, the p-value is
even larger (X2 = 0.113, p = 0.737). This is
also the case with RE: (X2 = 1.458, p =
0.2272) and (X2 = 0.098, p =
0.754).
Base Diameter
between Belowground Treatments for Sapindus
saponaria?Base Diameter
between Aboveground Treatments for Sapindus
saponaria?glmG.Diam.SASA2 <- glm(data = WetMassExp.SASA, Diameter.Base..cm. ~ Aboveground + Belowground,
family = Gamma(link = "log"))
glmG.Diam.C.SASA2 <- glm(data = WetMassExp.C.SASA, Diameter.Base..cm. ~ Aboveground + Belowground,
family = Gamma(link = "log"))
glmG.Diam.SASA2_glmmTMB <- glmmTMB(Diameter.Base..cm. ~ Aboveground + Belowground + (1 | Plot),
family = Gamma(link = "log"),
data = WetMassExp.SASA)
glmG.Diam.C.SASA2_glmmTMB <- glmmTMB(Diameter.Base..cm. ~ Aboveground + Belowground + (1 | Plot),
family = Gamma(link = "log"),
data = WetMassExp.C.SASA)
Anova(glmG.Diam.SASA2)
## Analysis of Deviance Table (Type II tests)
##
## Response: Diameter.Base..cm.
## LR Chisq Df Pr(>Chisq)
## Aboveground 12.761 1 0.000354 ***
## Belowground 3.178 1 0.074622 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(glmG.Diam.C.SASA2)
## Analysis of Deviance Table (Type II tests)
##
## Response: Diameter.Base..cm.
## LR Chisq Df Pr(>Chisq)
## Aboveground 21.253 1 4.03e-06 ***
## Belowground 0.933 1 0.334
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(glmG.Diam.SASA2_glmmTMB)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Diameter.Base..cm.
## Chisq Df Pr(>Chisq)
## Aboveground 32.501 1 1.19e-08 ***
## Belowground 4.629 1 0.0314 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(glmG.Diam.SASA2_glmmTMB)
## Family: Gamma ( log )
## Formula: Diameter.Base..cm. ~ Aboveground + Belowground + (1 | Plot)
## Data: WetMassExp.SASA
##
## AIC BIC logLik deviance df.resid
## 28.4 35.6 -9.2 18.4 26
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 0.181 0.426
## Number of obs: 31, groups: Plot, 6
##
## Dispersion estimate for Gamma family (sigma^2): 0.303
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.326 0.251 -5.28 1.3e-07 ***
## AbovegroundExposed 1.227 0.215 5.70 1.2e-08 ***
## BelowgroundUnprotected 0.461 0.214 2.15 0.031 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(glmG.Diam.C.SASA2_glmmTMB)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Diameter.Base..cm.
## Chisq Df Pr(>Chisq)
## Aboveground 43.891 1 3.47e-11 ***
## Belowground 0.481 1 0.488
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Answer:
Yes, there is a significant difference between
Aboveground Treatments on Base Diameter for
Sapindus saponaria (X2 = 12.761, p =
0.000354). When we exclude the outlier, this p-value is
even smaller (X2 = 21.253, p = 4.03e-06).
This is also the case with RE: (X2 = 32.501, p =
1.19e-08) and (X2 = 43.891, p =
3.47e-11).
Yes, there is a marginally significant difference between
Belowground Treatments on Base Diameter for
Sapindus saponaria (X2 = 3.178, p =
0.074622). When we exclude the outlier, we do not have
a difference (X2 = 0.933, p = 0.334). With
random effects, there is a statistically significant difference
(X2 = 4.629, p = 0.0314), but not when we
exclude the outlier (X2 = 0.481, p =
0.488).
Belowground Treatment
and Aboveground Treatment on Base Diameter for
Cymbopetalum baillonii?Anova(glmG.Diam.CYBA)
## Analysis of Deviance Table (Type II tests)
##
## Response: Diameter.Base..cm.
## LR Chisq Df Pr(>Chisq)
## Aboveground 6.136 1 0.0132 *
## Belowground 0.572 1 0.4494
## Aboveground:Belowground 1.694 1 0.1930
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(glmG.Diam.CYBA_glmmTMB)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Diameter.Base..cm.
## Chisq Df Pr(>Chisq)
## Aboveground 7.235 1 0.00715 **
## Belowground 0.737 1 0.39072
## Aboveground:Belowground 2.152 1 0.14240
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Answer: No, there is no interaction between
Belowground Treatment and
Aboveground Treatment on Base Diameter for
Cymbopetalum baillonii (X2 = 1.694, p =
0.193). Version with RE: (X2 = 2.152, p =
0.1424).
Base Diameter
between Belowground Treatments for Cymbopetalum
baillonii?Base Diameter
between Aboveground Treatments for Cymbopetalum
baillonii?glmG.Diam.CYBA2 <- glm(data = WetMassExp.CYBA, Diameter.Base..cm. ~ Aboveground + Belowground,
family = Gamma(link = "log"))
glmG.Diam.CYBA2_glmmTMB <- glmmTMB(Diameter.Base..cm. ~ Aboveground + Belowground + (1 | Plot),
family = Gamma(link = "log"),
data = WetMassExp.CYBA)
Anova(glmG.Diam.CYBA2)
## Analysis of Deviance Table (Type II tests)
##
## Response: Diameter.Base..cm.
## LR Chisq Df Pr(>Chisq)
## Aboveground 5.332 1 0.0209 *
## Belowground 0.497 1 0.4808
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(glmG.Diam.CYBA2_glmmTMB)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Diameter.Base..cm.
## Chisq Df Pr(>Chisq)
## Aboveground 6.696 1 0.00966 **
## Belowground 0.651 1 0.41962
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(glmG.Diam.CYBA2_glmmTMB)
## Family: Gamma ( log )
## Formula: Diameter.Base..cm. ~ Aboveground + Belowground + (1 | Plot)
## Data: WetMassExp.CYBA
##
## AIC BIC logLik deviance df.resid
## -14.1 -9.4 12.1 -24.1 14
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 7.76e-11 8.81e-06
## Number of obs: 19, groups: Plot, 5
##
## Dispersion estimate for Gamma family (sigma^2): 0.107
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.093 0.112 -9.79 <2e-16 ***
## AbovegroundExposed 0.423 0.163 2.59 0.0097 **
## BelowgroundUnprotected 0.123 0.152 0.81 0.4196
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Answer:
Yes, there is a significant difference between
Aboveground Treatments on Base Diameter for
Cymbopetalum baillonii (X2 = 5.332, p =
0.0209). Version with RE: (X2 = 6.696, p =
0.00966).
No, there is no significant difference between
Belowground Treatments on Base Diameter for
Cymbopetalum baillonii (X2 = 0.497, p =
0.4808). Version with RE: (X2 = 0.651, p =
0.41962).
Dry Mass between
species?glmDry_spp <- glm(data = DryMassExp, TotalDryMass ~ Species + Aboveground + Belowground,
family = Gamma(link = "log"))
glmDry_spp_glmmTMB <- glmmTMB(TotalDryMass ~ Species + Aboveground + Belowground + (1 | Plot),
family = Gamma(link = "log"),
data = DryMassExp)
Anova(glmDry_spp)
## Analysis of Deviance Table (Type II tests)
##
## Response: TotalDryMass
## LR Chisq Df Pr(>Chisq)
## Species 0.771 1 0.380
## Aboveground 27.073 1 1.96e-07 ***
## Belowground 2.255 1 0.133
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(glmDry_spp_glmmTMB)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: TotalDryMass
## Chisq Df Pr(>Chisq)
## Species 3.327 1 0.0681 .
## Aboveground 67.101 1 2.58e-16 ***
## Belowground 0.137 1 0.7110
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(glmDry_spp_glmmTMB)
## Family: Gamma ( log )
## Formula:
## TotalDryMass ~ Species + Aboveground + Belowground + (1 | Plot)
## Data: DryMassExp
##
## AIC BIC logLik deviance df.resid
## 209.1 220.5 -98.5 197.1 44
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 0.372 0.61
## Number of obs: 50, groups: Plot, 6
##
## Dispersion estimate for Gamma family (sigma^2): 0.902
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.632 0.429 -1.48 0.140
## SpeciesSapindus saponaria 0.572 0.314 1.82 0.068 .
## AbovegroundExposed 2.701 0.330 8.19 2.6e-16 ***
## BelowgroundProtected -0.118 0.319 -0.37 0.711
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Answer: No, there is not a significant difference
between Species on Dry Mass (X2 =
0.771, p = 0.380). Version with RE was marginally
significant (X2 = 3.327, p = 0.0681).
Belowground Treatment
and Aboveground Treatment on Dry Mass for
Sapindus saponaria?glmDry.SASA <- glm(data = DryMassExp.SASA, TotalDryMass ~ Aboveground * Belowground,
family = Gamma(link = "log"))
glmDry.SASA_glmmTMB <- glmmTMB(TotalDryMass ~ Aboveground * Belowground + (1 | Plot),
family = Gamma(link = "log"),
data = DryMassExp.SASA)
Anova(glmDry.SASA_glmmTMB)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: TotalDryMass
## Chisq Df Pr(>Chisq)
## Aboveground 51.611 1 6.77e-13 ***
## Belowground 1.779 1 0.182
## Aboveground:Belowground 0.399 1 0.528
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Answer: No, there was no significant interaction
between Belowground Treatment and
Aboveground Treatment on Dry Mass for
Sapindus saponaria (X2 = 0.495, p =
0.482). Version with RE: (X2 = 0.399, p =
0.528).
Dry Mass between
Belowground Treatments for Sapindus
saponaria?Dry Mass between
Aboveground Treatments for Sapindus
saponaria?glmDry.SASA2 <- glm(data = DryMassExp.SASA, TotalDryMass ~ Aboveground + Belowground,
family = Gamma(link = "log"))
glmDry.SASA2_glmmTMB <- glmmTMB(TotalDryMass ~ Aboveground + Belowground + (1 | Plot),
family = Gamma(link = "log"),
data = DryMassExp.SASA)
Anova(glmDry.SASA2)
## Analysis of Deviance Table (Type II tests)
##
## Response: TotalDryMass
## LR Chisq Df Pr(>Chisq)
## Aboveground 26.204 1 3.07e-07 ***
## Belowground 1.984 1 0.159
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(glmDry.SASA2_glmmTMB)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: TotalDryMass
## Chisq Df Pr(>Chisq)
## Aboveground 50.897 1 9.73e-13 ***
## Belowground 1.767 1 0.184
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(glmDry.SASA2_glmmTMB)
## Family: Gamma ( log )
## Formula: TotalDryMass ~ Aboveground + Belowground + (1 | Plot)
## Data: DryMassExp.SASA
##
## AIC BIC logLik deviance df.resid
## 152.7 159.9 -71.3 142.7 26
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 0.471 0.686
## Number of obs: 31, groups: Plot, 6
##
## Dispersion estimate for Gamma family (sigma^2): 0.891
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.0949 0.4547 0.21 0.83
## AbovegroundExposed 2.8559 0.4003 7.13 9.7e-13 ***
## BelowgroundProtected -0.5129 0.3859 -1.33 0.18
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Answer:
Yes, there is a significant difference between
Aboveground Treatments on Dry Mass for
Sapindus saponaria (X2 = 26.204, p =
3.07e-07). Version with RE: (X2 = 50.897, p
= 9.73e-13).
No, there was no significant difference between
Belowground Treatments on Dry Mass for
Sapindus saponaria (X2 = 1.984, p =
0.159). Version with RE: (X2 = 1.767, p =
0.184).
Belowground Treatment
and Aboveground Treatment on Dry Mass for
Cymbopetalum baillonii?glmDry.CYBA <- glm(data = DryMassExp.CYBA, TotalDryMass ~ Aboveground * Belowground,
family = Gamma(link = "log"))
glmDry.CYBA_glmmTMB <- glmmTMB(TotalDryMass ~ Aboveground * Belowground + (1 | Plot),
family = Gamma(link = "log"),
data = DryMassExp.CYBA)
Anova(glmDry.CYBA)
## Analysis of Deviance Table (Type II tests)
##
## Response: TotalDryMass
## LR Chisq Df Pr(>Chisq)
## Aboveground 12.953 1 0.000319 ***
## Belowground 0.732 1 0.392331
## Aboveground:Belowground 0.014 1 0.906206
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(glmDry.CYBA_glmmTMB)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: TotalDryMass
## Chisq Df Pr(>Chisq)
## Aboveground 12.533 1 0.0004 ***
## Belowground 0.780 1 0.3772
## Aboveground:Belowground 0.015 1 0.9031
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(glmDry.CYBA_glmmTMB)
## Family: Gamma ( log )
## Formula: TotalDryMass ~ Aboveground * Belowground + (1 | Plot)
## Data: DryMassExp.CYBA
##
## AIC BIC logLik deviance df.resid
## 59.2 64.9 -23.6 47.2 13
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 2.58e-10 1.61e-05
## Number of obs: 19, groups: Plot, 5
##
## Dispersion estimate for Gamma family (sigma^2): 0.704
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.00995 0.34261 0.03 0.977
## AbovegroundExposed 1.41556 0.59342 2.38 0.017 *
## BelowgroundProtected -0.37277 0.46690 -0.80 0.425
## AbovegroundExposed:BelowgroundProtected 0.10098 0.82917 0.12 0.903
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Answer: No, there was no significant interaction
among Belowground Treatment and
Aboveground Treatment on Dry Mass for
Cymbopetalum baillonii (X2 = 0.014, p =
0.906). Version with RE: (X2 = 0.015, p =
0.903)
Dry Mass between
Belowground Treatments for Cymbopetalum
baillonii?Dry Mass between
Aboveground Treatments for Cymbopetalum
baillonii?glmDry.CYBA2 <- glm(data = DryMassExp.CYBA, TotalDryMass ~ Aboveground + Belowground,
family = Gamma(link = "log"))
glmDry.CYBA2_glmmTMB <- glmmTMB(TotalDryMass ~ Aboveground + Belowground + (1 | Plot),
family = Gamma(link = "log"),
data = DryMassExp.CYBA)
Anova(glmDry.CYBA2)
## Analysis of Deviance Table (Type II tests)
##
## Response: TotalDryMass
## LR Chisq Df Pr(>Chisq)
## Aboveground 13.775 1 0.000206 ***
## Belowground 0.778 1 0.377709
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(glmDry.CYBA2_glmmTMB)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: TotalDryMass
## Chisq Df Pr(>Chisq)
## Aboveground 12.55 1 0.000397 ***
## Belowground 0.78 1 0.377002
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Answer:
Yes, there is a significant difference between
Aboveground Treatments on Dry Mass for
Cymbopetalum baillonii (X2 = 13.775, p =
0.0002). Version with RE: (X2 = 12.55, p =
0.0004).
No, there was no significant difference between
Belowground Treatments on Dry Mass for
Cymbopetalum baillonii (X2 = 0.778, p =
0.377709). Version with RE: (X2 = 0.78, p =
0.377).
Dry R-S Ratio between
species?lm5Dry_spp <- lm(data = DryMassExp, Ratio ~ Species + Aboveground + Belowground)
lm5Dry_spp_RE <- lmer(data = DryMassExp, Ratio ~ Species + Aboveground + Belowground + (1 | Plot))
Anova(lm5Dry_spp)
## Anova Table (Type II tests)
##
## Response: Ratio
## Sum Sq Df F value Pr(>F)
## Species 1.651 1 9.477 0.0035 **
## Aboveground 0.627 1 3.599 0.0641 .
## Belowground 0.586 1 3.367 0.0730 .
## Residuals 8.013 46
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(lm5Dry_spp_RE)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Ratio
## Chisq Df Pr(>Chisq)
## Species 9.477 1 0.00208 **
## Aboveground 3.599 1 0.05780 .
## Belowground 3.367 1 0.06652 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(lm5Dry_spp_RE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Ratio ~ Species + Aboveground + Belowground + (1 | Plot)
## Data: DryMassExp
##
## REML criterion at convergence: 61.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.445 -0.596 -0.159 0.437 3.347
##
## Random effects:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 0.000 0.000
## Residual 0.174 0.417
## Number of obs: 50, groups: Plot, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.504 0.120 4.19
## SpeciesSapindus saponaria 0.382 0.124 3.08
## AbovegroundExposed 0.230 0.121 1.90
## BelowgroundProtected 0.217 0.118 1.83
##
## Correlation of Fixed Effects:
## (Intr) SpcsSs AbvgrE
## SpcsSpndssp -0.541
## AbvgrndExps -0.315 -0.196
## BlwgrndPrtc -0.516 -0.020 -0.006
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
Answer: Yes, there is a significant difference
between Species on Dry R-S Ratios
(F1,46 = 9.477, p = 0.0035). Version with
RE: (F1,46 = 9.477, p = 0.00208).
Belowground Treatment
and Aboveground Treatment on
Dry Root-to-Shoot Ratios for Sapindus
saponaria?lm5Dry.SASA <- lm(data = DryMassExp.SASA, Ratio ~ Aboveground * Belowground)
lm5Dry.SASA_RE <- lmer(data = DryMassExp, Ratio ~ Species + Aboveground + Belowground + (1 | Plot))
Anova(lm5Dry.SASA)
## Anova Table (Type II tests)
##
## Response: Ratio
## Sum Sq Df F value Pr(>F)
## Aboveground 0.407 1 1.878 0.182
## Belowground 0.453 1 2.092 0.160
## Aboveground:Belowground 0.003 1 0.014 0.908
## Residuals 5.851 27
Anova(lm5Dry.SASA_RE)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Ratio
## Chisq Df Pr(>Chisq)
## Species 9.477 1 0.00208 **
## Aboveground 3.599 1 0.05780 .
## Belowground 3.367 1 0.06652 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Answer: No, there was no significant interaction
between Belowground Treatment and
Aboveground Treatment on Dry R-S Ratios for
Sapindus saponaria (F1,27 = 0.014, p = 0.908).
Version with RE showed a marginally significant interaction:
(F1,27 = 3.367, p = 0.067).
Dry R-S Ratios
between Belowground Treatments for Sapindus
saponaria?Dry R-S Ratios
between Aboveground Treatments for Sapindus
saponaria?lm5Dry.SASA2 <- lm(data = DryMassExp.SASA, Ratio ~ Aboveground + Belowground)
lm5Dry.SASA2_RE <- lmer(data = DryMassExp.SASA, Ratio ~ Aboveground + Belowground + (1 | Plot))
Anova(lm5Dry.SASA2)
## Anova Table (Type II tests)
##
## Response: Ratio
## Sum Sq Df F value Pr(>F)
## Aboveground 0.407 1 1.947 0.174
## Belowground 0.453 1 2.169 0.152
## Residuals 5.854 28
Anova(lm5Dry.SASA2_RE)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Ratio
## Chisq Df Pr(>Chisq)
## Aboveground 1.947 1 0.163
## Belowground 2.169 1 0.141
summary(lm5Dry.SASA2_RE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Ratio ~ Aboveground + Belowground + (1 | Plot)
## Data: DryMassExp.SASA
##
## REML criterion at convergence: 43.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.345 -0.680 -0.135 0.466 3.029
##
## Random effects:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 0.000 0.000
## Residual 0.209 0.457
## Number of obs: 31, groups: Plot, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.872 0.147 5.92
## AbovegroundExposed 0.229 0.164 1.40
## BelowgroundProtected 0.243 0.165 1.47
##
## Correlation of Fixed Effects:
## (Intr) AbvgrE
## AbvgrndExps -0.558
## BlwgrndPrtc -0.598 -0.029
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
Answer:
No, there is no significant difference between
Aboveground Treatments on Dry R-S Ratios for
Sapindus saponaria (F1,28 = 2.169, p =
0.174). Version with RE: (F1,28 = 1.947, p =
0.163).
No, there is no marginally significant difference between
Belowground Treatments on Dry R-S Ratios for
Sapindus saponaria (F1,28 = 2.169, p =
0.152). Version with RE: (F1,28 = 2.169, p =
0.141).
Belowground Treatment
and Aboveground Treatment on
Dry Root-to-Shoot Ratios for Cymbopetalum
baillonii?lm5Dry.CYBA <- lm(data = DryMassExp.CYBA, Ratio ~ Aboveground * Belowground)
lm5Dry.CYBA_RE <- lmer(data = DryMassExp.CYBA, Ratio ~ Aboveground * Belowground + (1 | Plot))
Anova(lm5Dry.CYBA)
## Anova Table (Type II tests)
##
## Response: Ratio
## Sum Sq Df F value Pr(>F)
## Aboveground 0.2136 1 1.645 0.219
## Belowground 0.1458 1 1.123 0.306
## Aboveground:Belowground 0.1974 1 1.520 0.237
## Residuals 1.9481 15
Anova(lm5Dry.CYBA_RE)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Ratio
## Chisq Df Pr(>Chisq)
## Aboveground 1.390 1 0.238
## Belowground 0.792 1 0.373
## Aboveground:Belowground 1.610 1 0.204
Answer: No, there was no significant interaction
between Belowground Treatment and
Aboveground Treatment on Dry R-S Ratios for
Cymbopetalum baillonii (F1,15 = 1.520, p =
0.237). Version with RE: (F1,15 = 1.610, p =
0.204).
Dry R-S Ratios
between Belowground Treatments for Cymbopetalum
baillonii?Dry R-S Ratios
between Aboveground Treatments for Cymbopetalum
baillonii?lm5Dry.CYBA2 <- lm(data = DryMassExp.CYBA, Ratio ~ Aboveground + Belowground)
lm5Dry.CYBA2_RE <- lmer(data = DryMassExp.CYBA, Ratio ~ Aboveground + Belowground + (1 | Plot))
Anova(lm5Dry.CYBA2)
## Anova Table (Type II tests)
##
## Response: Ratio
## Sum Sq Df F value Pr(>F)
## Aboveground 0.2136 1 1.593 0.225
## Belowground 0.1458 1 1.087 0.313
## Residuals 2.1455 16
Anova(lm5Dry.CYBA2_RE)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Ratio
## Chisq Df Pr(>Chisq)
## Aboveground 1.350 1 0.245
## Belowground 0.787 1 0.375
summary(lm5Dry.CYBA2_RE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Ratio ~ Aboveground + Belowground + (1 | Plot)
## Data: DryMassExp.CYBA
##
## REML criterion at convergence: 18.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.271 -0.604 -0.154 0.402 2.104
##
## Random effects:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 0.025 0.158
## Residual 0.117 0.343
## Number of obs: 19, groups: Plot, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.536 0.164 3.26
## AbovegroundExposed 0.245 0.211 1.16
## BelowgroundProtected 0.150 0.169 0.89
##
## Correlation of Fixed Effects:
## (Intr) AbvgrE
## AbvgrndExps -0.523
## BlwgrndPrtc -0.577 0.065
Answer:
No, there is no significant difference between
Aboveground Treatments on Dry R-S Ratios for
Cymbopetalum baillonii (F1,16 = 1.593, p =
0.225). Version with RE: (F1,16 = 1.350, p =
0.245).
No, there is no marginally significant difference between
Belowground Treatments on Dry R-S Ratios for
Cymbopetalum baillonii (F1,16 = 1.087, p =
0.313). Version with RE: (F1,16 = 0.787, p =
0.375).
Leaf Dry Mass between
Species?glm.nb.Leaf_spp <- glm(data = DryMassExp, LeafAveMass ~ Species + Aboveground + Belowground,
family = Gamma(link = "log"))
glm.nb.Leaf_spp_RE <- glmmTMB(LeafAveMass ~ Species + Aboveground + Belowground + (1 | Plot),
family = Gamma(link = "log"),
data = DryMassExp)
Anova(glm.nb.Leaf_spp)
## Analysis of Deviance Table (Type II tests)
##
## Response: LeafAveMass
## LR Chisq Df Pr(>Chisq)
## Species 2.82 1 0.0932 .
## Aboveground 47.10 1 6.75e-12 ***
## Belowground 1.28 1 0.2572
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(glm.nb.Leaf_spp_RE)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: LeafAveMass
## Chisq Df Pr(>Chisq)
## Species 4.029 1 0.0447 *
## Aboveground 52.051 1 5.41e-13 ***
## Belowground 0.421 1 0.5165
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(glm.nb.Leaf_spp_RE)
## Family: Gamma ( log )
## Formula:
## LeafAveMass ~ Species + Aboveground + Belowground + (1 | Plot)
## Data: DryMassExp
##
## AIC BIC logLik deviance df.resid
## -165.0 -153.6 88.5 -177.0 44
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 0.0657 0.256
## Number of obs: 50, groups: Plot, 6
##
## Dispersion estimate for Gamma family (sigma^2): 0.469
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.608 0.255 -14.15 < 2e-16 ***
## SpeciesSapindus saponaria 0.438 0.218 2.01 0.045 *
## AbovegroundExposed 1.641 0.228 7.21 5.4e-13 ***
## BelowgroundProtected -0.136 0.209 -0.65 0.516
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Answer: There is a marginally significant difference
between Species on Leaf Dry Mass
(F1,46 = 2.82, p = 0.0932). With RE there is
a marginally significant difference (F1,46 = 2.955, p =
0.0856).
Belowground Treatment
and Aboveground Treatment on Leaf Dry Mass for
Sapindus saponaria?glm.nb.Leaf.SASA <- glm(data = DryMassExp.SASA, LeafAveMass ~ Aboveground * Belowground,
family = Gamma(link = "log"))
glm.nb.Leaf.SASA_RE <- glmmTMB(LeafAveMass ~ Aboveground * Belowground + (1 | Plot),
family = Gamma(link = "log"),
data = DryMassExp.SASA)
Anova(glm.nb.Leaf.SASA)
## Analysis of Deviance Table (Type II tests)
##
## Response: LeafAveMass
## LR Chisq Df Pr(>Chisq)
## Aboveground 38.18 1 6.46e-10 ***
## Belowground 1.04 1 0.309
## Aboveground:Belowground 0.07 1 0.796
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(glm.nb.Leaf.SASA_RE)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: LeafAveMass
## Chisq Df Pr(>Chisq)
## Aboveground 58.770 1 1.77e-14 ***
## Belowground 1.136 1 0.287
## Aboveground:Belowground 0.100 1 0.752
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Answer: No, there was no significant interaction
between Belowground Treatment and
Aboveground Treatment on Leaf Dry Mass for
Sapindus saponaria (X2 = 0.07, p = 0.796).
Version with RE: (X2 = 0.100, p = 0.752).
Leaf Dry Mass
between Aboveground Treatments for Sapindus
saponaria?Leaf Dry Mass
between Belowground Treatments for Sapindus
saponaria?glm.nb.Leaf.SASA2 <- glm(data = DryMassExp.SASA, LeafAveMass ~ Aboveground + Belowground,
family = Gamma(link = "log"))
glm.nb.Leaf.SASA2_RE <- glmmTMB(LeafAveMass ~ Aboveground + Belowground + (1 | Plot),
family = Gamma(link = "log"),
data = DryMassExp.SASA)
Anova(glm.nb.Leaf.SASA2)
## Analysis of Deviance Table (Type II tests)
##
## Response: LeafAveMass
## LR Chisq Df Pr(>Chisq)
## Aboveground 39.12 1 3.98e-10 ***
## Belowground 1.06 1 0.303
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(glm.nb.Leaf.SASA2_RE)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: LeafAveMass
## Chisq Df Pr(>Chisq)
## Aboveground 58.709 1 1.83e-14 ***
## Belowground 1.135 1 0.287
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(glm.nb.Leaf.SASA2_RE)
## Family: Gamma ( log )
## Formula: LeafAveMass ~ Aboveground + Belowground + (1 | Plot)
## Data: DryMassExp.SASA
##
## AIC BIC logLik deviance df.resid
## -86.0 -78.9 48.0 -96.0 26
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 0.0847 0.291
## Number of obs: 31, groups: Plot, 6
##
## Dispersion estimate for Gamma family (sigma^2): 0.399
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.207 0.251 -12.79 < 2e-16 ***
## AbovegroundExposed 1.830 0.239 7.66 1.8e-14 ***
## BelowgroundProtected -0.251 0.236 -1.07 0.29
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Answer:
Yes, there is a significant difference between
Aboveground Treatments on Leaf Dry Mass for
Sapindus saponaria (X2 = 39.12, p = 3.98e-10).
Version with RE: (X2 = 58.709, p = 1.83e-14).
No, there is no significant difference between
Belowground Treatments on Leaf Dry Mass for
Sapindus saponaria (X2 = 1.06, p = 0.303).
Version with RE: (F1,28 = 1.135, p =
0.287).
Belowground Treatment
and Aboveground Treatment on Leaf Dry Mass for
Cymbopetalum baillonii?glm.nb.Leaf.CYBA <- glm(data = DryMassExp.CYBA, LeafAveMass ~ Aboveground * Belowground,
family = Gamma(link = "log"))
glm.nb.Leaf.CYBA_RE <- glmmTMB(LeafAveMass ~ Aboveground * Belowground + (1 | Plot),
family = Gamma(link = "log"),
data = DryMassExp.CYBA)
Anova(glm.nb.Leaf.CYBA)
## Analysis of Deviance Table (Type II tests)
##
## Response: LeafAveMass
## LR Chisq Df Pr(>Chisq)
## Aboveground 8.821 1 0.00298 **
## Belowground 0.389 1 0.53286
## Aboveground:Belowground 0.228 1 0.63326
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(glm.nb.Leaf.CYBA_RE)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: LeafAveMass
## Chisq Df Pr(>Chisq)
## Aboveground 8.927 1 0.00281 **
## Belowground 0.433 1 0.51049
## Aboveground:Belowground 0.252 1 0.61601
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Answer: No, there was no significant interaction
between Belowground Treatment and
Aboveground Treatment on Leaf Dry Mass for
Sapindus saponaria (X2 = 0.228, p = 0.63326).
Version with RE: (X2 = 0.252, p = 0.616).
Leaf Dry Mass
between Aboveground Treatments for Cymbopetalum
baillonii?Leaf Dry Mass
between Belowground Treatments for Cymbopetalum
baillonii?glm.nb.Leaf.CYBA2 <- glm(data = DryMassExp.CYBA, LeafAveMass ~ Aboveground + Belowground,
family = Gamma(link = "log"))
glm.nb.Leaf.CYBA2_RE <- glmmTMB(LeafAveMass ~ Aboveground + Belowground + (1 | Plot),
family = Gamma(link = "log"),
data = DryMassExp.CYBA)
Anova(glm.nb.Leaf.CYBA2)
## Analysis of Deviance Table (Type II tests)
##
## Response: LeafAveMass
## LR Chisq Df Pr(>Chisq)
## Aboveground 8.825 1 0.00297 **
## Belowground 0.389 1 0.53277
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(glm.nb.Leaf.CYBA2_RE)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: LeafAveMass
## Chisq Df Pr(>Chisq)
## Aboveground 8.810 1 0.003 **
## Belowground 0.423 1 0.516
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(glm.nb.Leaf.CYBA2_RE)
## Family: Gamma ( log )
## Formula: LeafAveMass ~ Aboveground + Belowground + (1 | Plot)
## Data: DryMassExp.CYBA
##
## AIC BIC logLik deviance df.resid
## -74.1 -69.4 42.1 -84.1 14
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 2.86e-10 1.69e-05
## Number of obs: 19, groups: Plot, 5
##
## Dispersion estimate for Gamma family (sigma^2): 0.531
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.338 0.278 -12.01 <2e-16 ***
## AbovegroundExposed 1.075 0.362 2.97 0.003 **
## BelowgroundProtected -0.219 0.337 -0.65 0.516
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Answer:
Yes, there is a significant difference between
Aboveground Treatments on Base Diameter for
Cymbopetalum baillonii (F1,16 = 8.825, p =
0.00297). Version with RE was statistically
significant: (X2 = 8.810, p = 0.003).
No, there is no significant difference between
Belowground Treatments on Base Diameter for
Cymbopetalum baillonii (x2 = 0.389, p =
0.533). Version with RE (X2 = 0.423, p =
0.516).
MoistureVOrig <- lmer(Moisture_Leaves2 ~ Species * Aboveground * Belowground + (1 | Plot), data = DryMassExp)
Anova(MoistureVOrig)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Moisture_Leaves2
## Chisq Df Pr(>Chisq)
## Species 11.781 1 0.000598 ***
## Aboveground 0.488 1 0.484958
## Belowground 1.286 1 0.256753
## Species:Aboveground 0.127 1 0.721596
## Species:Belowground 0.628 1 0.428250
## Aboveground:Belowground 3.527 1 0.060394 .
## Species:Aboveground:Belowground 1.415 1 0.234270
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
There is no significant interaction among the Species, Aboveground, or Belowground Treatments (X2 = 1.415, p = 0.234).
Moisture Content (%) between
Species?I want to check the model diagnostics for a normal model.
plot(lm(data = DryMassExp, Moisture2 ~ Aboveground * Belowground + Species))
As I thought, the normal Gaussian model does not behave in odd ways. It
might be an overkill to use beta regression after all.
betaM.1_spp_glmmTMB <- glmmTMB(Moisture2 ~ Species + Aboveground + Belowground + (1 | Plot),
family = beta_family(),
data = DryMassExp)
Anova(betaM.1_spp_glmmTMB)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Moisture2
## Chisq Df Pr(>Chisq)
## Species 39.091 1 4.05e-10 ***
## Aboveground 4.184 1 0.0408 *
## Belowground 0.008 1 0.9269
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(betaM.1_spp_glmmTMB)
## Family: beta ( logit )
## Formula: Moisture2 ~ Species + Aboveground + Belowground + (1 | Plot)
## Data: DryMassExp
##
## AIC BIC logLik deviance df.resid
## -84.3 -72.8 48.1 -96.3 44
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 0.0153 0.124
## Number of obs: 50, groups: Plot, 6
##
## Dispersion parameter for beta family (): 28.5
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.7571 0.1401 -5.40 6.6e-08 ***
## SpeciesSapindus saponaria 0.7182 0.1149 6.25 4.0e-10 ***
## AbovegroundExposed 0.2731 0.1335 2.05 0.041 *
## BelowgroundProtected 0.0107 0.1170 0.09 0.927
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
betaM.1V2_spp_glmmTMB <- glmmTMB(Moisture_NoRoots2 ~ Species + Aboveground + Belowground + (1 | Plot),
family = beta_family(),
data = DryMassExp)
Anova(betaM.1V2_spp_glmmTMB)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Moisture_NoRoots2
## Chisq Df Pr(>Chisq)
## Species 34.725 1 3.8e-09 ***
## Aboveground 4.351 1 0.037 *
## Belowground 0.135 1 0.713
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(betaM.1V2_spp_glmmTMB)
## Family: beta ( logit )
## Formula:
## Moisture_NoRoots2 ~ Species + Aboveground + Belowground + (1 | Plot)
## Data: DryMassExp
##
## AIC BIC logLik deviance df.resid
## -80.7 -69.2 46.4 -92.7 44
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 0.0126 0.112
## Number of obs: 50, groups: Plot, 6
##
## Dispersion parameter for beta family (): 25.6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.8671 0.1468 -5.91 3.5e-09 ***
## SpeciesSapindus saponaria 0.7203 0.1222 5.89 3.8e-09 ***
## AbovegroundExposed 0.2978 0.1427 2.09 0.037 *
## BelowgroundProtected -0.0453 0.1231 -0.37 0.713
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
betaM.1V3_spp_glmmTMB <- glmmTMB(Moisture_Leaves2 ~ Species + Aboveground + Belowground + (1 | Plot),
family = beta_family(),
data = DryMassExp)
Anova(betaM.1V3_spp_glmmTMB)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Moisture_Leaves2
## Chisq Df Pr(>Chisq)
## Species 11.452 1 0.000714 ***
## Aboveground 0.995 1 0.318419
## Belowground 0.683 1 0.408592
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(betaM.1V3_spp_glmmTMB)
## Family: beta ( logit )
## Formula:
## Moisture_Leaves2 ~ Species + Aboveground + Belowground + (1 | Plot)
## Data: DryMassExp
##
## AIC BIC logLik deviance df.resid
## -54.0 -42.5 33.0 -66.0 44
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 0.0506 0.225
## Number of obs: 50, groups: Plot, 6
##
## Dispersion parameter for beta family (): 14.5
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.834 0.203 -4.10 4.2e-05 ***
## SpeciesSapindus saponaria 0.554 0.164 3.38 0.00071 ***
## AbovegroundExposed 0.174 0.174 1.00 0.31842
## BelowgroundProtected -0.138 0.167 -0.83 0.40859
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Answer: Yes, there is a significant difference
between Species on Moisture Content (%)
(X2 = 37.157, p = 1.09e-09). Version with
RE: (X2 = 39.091, p = 4.05e-10). The same is
the case in the version without root mass. The same is the case with
only leaves. Let’s check a normal model for only leaves since this is
the one we will likely be using.
MoistureV3 <- lmer(Moisture_Leaves2 ~ Species +
Aboveground + Belowground + (1 | Plot), data = DryMassExp)
Anova(MoistureV3)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Moisture_Leaves2
## Chisq Df Pr(>Chisq)
## Species 10.715 1 0.00106 **
## Aboveground 0.534 1 0.46489
## Belowground 1.105 1 0.29309
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(MoistureV3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Moisture_Leaves2 ~ Species + Aboveground + Belowground + (1 |
## Plot)
## Data: DryMassExp
##
## REML criterion at convergence: -53.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.979 -0.609 -0.140 0.409 3.618
##
## Random effects:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 0.0034 0.0583
## Residual 0.0127 0.1125
## Number of obs: 50, groups: Plot, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.3212 0.0429 7.50
## SpeciesSapindus saponaria 0.1121 0.0343 3.27
## AbovegroundExposed 0.0264 0.0361 0.73
## BelowgroundProtected -0.0356 0.0338 -1.05
##
## Correlation of Fixed Effects:
## (Intr) SpcsSs AbvgrE
## SpcsSpndssp -0.433
## AbvgrndExps -0.317 -0.192
## BlwgrndPrtc -0.444 0.001 0.050
Answer: Yes, according to the normal model there is
a significant difference between Species on
Leaf Moisture Content (%) (X2 = 10.7146, p =
0.001063).
Belowground Treatment
and Aboveground Treatment on
Moisture Content (%) for Sapindus saponaria?betaM.1.SASA <- betareg(data = DryMassExp.SASA, Moisture2 ~ Aboveground * Belowground)
betaM.1.SASA_glmmTMB <- glmmTMB(Moisture2 ~ Aboveground * Belowground + (1 | Plot),
family = beta_family(),
data = DryMassExp.SASA)
Anova(betaM.1.SASA_glmmTMB)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Moisture2
## Chisq Df Pr(>Chisq)
## Aboveground 4.239 1 0.0395 *
## Belowground 2.343 1 0.1259
## Aboveground:Belowground 1.644 1 0.1997
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
betaM.1V2.SASA <- betareg(data = DryMassExp.SASA, Moisture_NoRoots2 ~ Aboveground * Belowground)
betaM.1V2.SASA_glmmTMB <- glmmTMB(Moisture_NoRoots2 ~ Aboveground * Belowground + (1 | Plot),
family = beta_family(),
data = DryMassExp.SASA)
Anova(betaM.1V2.SASA_glmmTMB)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Moisture_NoRoots2
## Chisq Df Pr(>Chisq)
## Aboveground 6.081 1 0.0137 *
## Belowground 4.445 1 0.0350 *
## Aboveground:Belowground 1.678 1 0.1952
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
betaM.1V2.SASA <- betareg(data = DryMassExp.SASA, Moisture_Leaves2 ~ Aboveground * Belowground)
betaM.1V2.SASA_glmmTMB <- glmmTMB(Moisture_Leaves2 ~ Aboveground * Belowground + (1 | Plot),
family = beta_family(),
data = DryMassExp.SASA)
Anova(betaM.1V2.SASA_glmmTMB)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Moisture_Leaves2
## Chisq Df Pr(>Chisq)
## Aboveground 1.154 1 0.283
## Belowground 1.916 1 0.166
## Aboveground:Belowground 1.261 1 0.261
Answer: No, there was no significant interaction
between Belowground Treatment and
Aboveground Treatment on Moisture Content for
Sapindus saponaria (X2 = 1.058, p =
0.304). Version with RE: (X2 = 1.644, p =
0.1997). Same goes for the model without root mass and
only leaves.
MoistureV3.2 <- lmer(Moisture_Leaves2 ~ Species +
Aboveground * Belowground + (1 | Plot), data = DryMassExp)
Anova(MoistureV3)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Moisture_Leaves2
## Chisq Df Pr(>Chisq)
## Species 10.715 1 0.00106 **
## Aboveground 0.534 1 0.46489
## Belowground 1.105 1 0.29309
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Answer: No, according to the normal model there is
no significant interaction Belowground Treatment and
Aboveground Treatment on
Leaf Moisture Content (%) (X2 = 1.1054, p =
0.293).
Moisture Content
between Aboveground Treatments for Sapindus
saponaria?#betaM.1.SASA2 <- betareg(data = DryMassExp.SASA, Moisture2 ~ Aboveground + Belowground)
betaM.1.SASA2_glmmTMB <- glmmTMB(Moisture2 ~ Aboveground + Belowground + (1 | Plot),
family = beta_family(),
data = DryMassExp.SASA)
#Anova(betaM.1.SASA2)
Anova(betaM.1.SASA2_glmmTMB)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Moisture2
## Chisq Df Pr(>Chisq)
## Aboveground 4.051 1 0.0441 *
## Belowground 2.147 1 0.1428
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(betaM.1.SASA2_glmmTMB)
## Family: beta ( logit )
## Formula: Moisture2 ~ Aboveground + Belowground + (1 | Plot)
## Data: DryMassExp.SASA
##
## AIC BIC logLik deviance df.resid
## -59.7 -52.6 34.9 -69.7 26
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 0.0214 0.146
## Number of obs: 31, groups: Plot, 6
##
## Dispersion parameter for beta family (): 45.6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.0753 0.1183 0.64 0.524
## AbovegroundExposed 0.2223 0.1104 2.01 0.044 *
## BelowgroundProtected -0.1659 0.1132 -1.47 0.143
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#betaM.1V2.SASA2 <- betareg(data = DryMassExp.SASA, Moisture_NoRoots2 ~ Aboveground + Belowground)
betaM.1V2.SASA2_glmmTMB <- glmmTMB(Moisture_NoRoots2 ~ Aboveground + Belowground + (1 | Plot),
family = beta_family(),
data = DryMassExp.SASA)
#Anova(betaM.1V2.SASA2)
Anova(betaM.1V2.SASA2_glmmTMB)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Moisture_NoRoots2
## Chisq Df Pr(>Chisq)
## Aboveground 5.783 1 0.0162 *
## Belowground 4.083 1 0.0433 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(betaM.1V2.SASA2_glmmTMB)
## Family: beta ( logit )
## Formula: Moisture_NoRoots2 ~ Aboveground + Belowground + (1 | Plot)
## Data: DryMassExp.SASA
##
## AIC BIC logLik deviance df.resid
## -56.5 -49.3 33.2 -66.5 26
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 0.0255 0.16
## Number of obs: 31, groups: Plot, 6
##
## Dispersion parameter for beta family (): 41.1
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.0445 0.1251 -0.36 0.722
## AbovegroundExposed 0.2799 0.1164 2.40 0.016 *
## BelowgroundProtected -0.2385 0.1180 -2.02 0.043 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
options(digits = 7)
#betaM.1V3.SASA2 <- betareg(data = DryMassExp.SASA, Moisture_Leaves2 ~ Aboveground + Belowground)
betaM.1V3.SASA2_glmmTMB <- glmmTMB(Moisture_Leaves2 ~ Aboveground + Belowground + (1 | Plot),
family = beta_family(),
data = DryMassExp.SASA)
#Anova(betaM.1V3.SASA2)
Anova(betaM.1V3.SASA2_glmmTMB)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Moisture_Leaves2
## Chisq Df Pr(>Chisq)
## Aboveground 1.1425 1 0.2851
## Belowground 1.7947 1 0.1804
summary(betaM.1V3.SASA2_glmmTMB)
## Family: beta ( logit )
## Formula: Moisture_Leaves2 ~ Aboveground + Belowground + (1 | Plot)
## Data: DryMassExp.SASA
##
## AIC BIC logLik deviance df.resid
## -29.9 -22.7 20.0 -39.9 26
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 0.03226 0.1796
## Number of obs: 31, groups: Plot, 6
##
## Dispersion parameter for beta family (): 15.3
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.2275 0.1877 -1.212 0.225
## AbovegroundExposed 0.1988 0.1860 1.069 0.285
## BelowgroundProtected -0.2626 0.1960 -1.340 0.180
Answer:
Yes, there is a significant difference between
Aboveground Treatments on Moisture Content for
Sapindus saponaria (X2 = 4.474, p =
0.0344). Version with RE: (X2 = 4.051, p =
0.0441). The version with no roots (X2 =
5.7833, p = 0.01618). The version with only leaves
(X2 = 1.7947, p = 0.1804).
No, there is no significant difference between
Belowground Treatments on Moisture Content for
Sapindus saponaria (X2 = 0.691, p =
0.4058). Version with RE: (X2 = 2.147, p =
0.1428), but the version without root mass did show a
significant difference (X2 = 4.083, p =
0.0433). The version with only leaves (X2 =
0.9107, p = 0.2838).
MoistureV3.3 <- lmer(Moisture_Leaves2 ~ Aboveground + Belowground + (1 | Plot), data = DryMassExp.SASA)
Anova(MoistureV3.3)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Moisture_Leaves2
## Chisq Df Pr(>Chisq)
## Aboveground 0.7927 1 0.3733
## Belowground 1.5093 1 0.2192
summary(MoistureV3.3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Moisture_Leaves2 ~ Aboveground + Belowground + (1 | Plot)
## Data: DryMassExp.SASA
##
## REML criterion at convergence: -30.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.3075 -0.5346 -0.2240 0.1771 3.4557
##
## Random effects:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 0.002305 0.04802
## Residual 0.013680 0.11696
## Number of obs: 31, groups: Plot, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.43662 0.04386 9.955
## AbovegroundExposed 0.03843 0.04317 0.890
## BelowgroundProtected -0.05345 0.04351 -1.229
##
## Correlation of Fixed Effects:
## (Intr) AbvgrE
## AbvgrndExps -0.509
## BlwgrndPrtc -0.539 -0.004
Answer: No, according to the normal model there is
no significant difference between Aboveground Treatment on
Leaf Moisture Content (%) (X2 = 0.7927, p =
0.3733).
Answer: No, according to the normal model there is
no significant difference between Belowground Treatment on
Leaf Moisture Content (%) (X2 = 1.5093, p =
0.2192).
Belowground Treatment
and Aboveground Treatment on
Moisture Content (%) for Cymbopetalum
baillonii?betaM.1.CYBA <- betareg(data = DryMassExp.CYBA, Moisture2 ~ Aboveground * Belowground)
betaM.1.CYBA_glmmTMB <- glmmTMB(Moisture2 ~ Aboveground * Belowground + (1 | Plot),
family = beta_family(),
data = DryMassExp.CYBA)
Anova(betaM.1.CYBA)
## Analysis of Deviance Table (Type II tests)
##
## Response: Moisture2
## Df Chisq Pr(>Chisq)
## Aboveground 1 8.3005 0.003963 **
## Belowground 1 3.2182 0.072825 .
## Aboveground:Belowground 1 2.9263 0.087149 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(betaM.1.CYBA_glmmTMB)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Moisture2
## Chisq Df Pr(>Chisq)
## Aboveground 8.3005 1 0.003963 **
## Belowground 3.2182 1 0.072825 .
## Aboveground:Belowground 2.9263 1 0.087149 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
betaM.1V2.CYBA <- betareg(data = DryMassExp.CYBA, Moisture_NoRoots2 ~ Aboveground * Belowground)
betaM.1V2.CYBA_glmmTMB <- glmmTMB(Moisture_NoRoots2 ~ Aboveground * Belowground + (1 | Plot),
family = beta_family(),
data = DryMassExp.CYBA)
Anova(betaM.1V2.CYBA)
## Analysis of Deviance Table (Type II tests)
##
## Response: Moisture_NoRoots2
## Df Chisq Pr(>Chisq)
## Aboveground 1 5.8581 0.01551 *
## Belowground 1 2.2468 0.13389
## Aboveground:Belowground 1 5.3695 0.02049 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(betaM.1V2.CYBA_glmmTMB)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Moisture_NoRoots2
## Chisq Df Pr(>Chisq)
## Aboveground 5.8581 1 0.01551 *
## Belowground 2.2468 1 0.13389
## Aboveground:Belowground 5.3695 1 0.02049 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
betaM.1V3.CYBA <- betareg(data = DryMassExp.CYBA, Moisture_Leaves2 ~ Aboveground * Belowground)
betaM.1V3.CYBA_glmmTMB <- glmmTMB(Moisture_Leaves2 ~ Aboveground * Belowground + (1 | Plot),
family = beta_family(),
data = DryMassExp.CYBA)
Anova(betaM.1V3.CYBA)
## Analysis of Deviance Table (Type II tests)
##
## Response: Moisture_Leaves2
## Df Chisq Pr(>Chisq)
## Aboveground 1 3.6035 0.05766 .
## Belowground 1 1.3247 0.24975
## Aboveground:Belowground 1 4.9001 0.02686 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(betaM.1V3.CYBA_glmmTMB)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Moisture_Leaves2
## Chisq Df Pr(>Chisq)
## Aboveground 3.6035 1 0.05766 .
## Belowground 1.3247 1 0.24975
## Aboveground:Belowground 4.9001 1 0.02686 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Answer: No, there was no significant interaction
between Belowground Treatment and
Aboveground Treatment on Moisture Content for
Cymbopetalum baillonii (X2 = 2.926, p =
0.088). Version with RE: (X2 = 2.926, p =
0.087). But the version without root mass did show an
interaction (X2 = 5.370, p = 0.0205), as did
the version with onyl leaves.
MoistureV3.4 <- lmer(Moisture_Leaves ~ Aboveground * Belowground + (1 | Plot), data = DryMassExp.CYBA)
Anova(MoistureV3.4)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Moisture_Leaves
## Chisq Df Pr(>Chisq)
## Aboveground 0.1057 1 0.745072
## Belowground 0.0749 1 0.784379
## Aboveground:Belowground 8.1004 1 0.004426 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Answer: There is also a significant interaction in the normal model (X2 = 8.1004, p = 0.004426).
Moisture Content
between Aboveground Treatments for Cymbopetalum
baillonii?Moisture Content
between Belowground Treatments for Cymbopetalum
baillonii?betaM.1.CYBA2 <- betareg(data = DryMassExp.CYBA, Moisture2 ~ Aboveground + Belowground)
betaM.1.CYBA2_glmmTMB <- glmmTMB(Moisture2 ~ Aboveground + Belowground + (1 | Plot),
family = beta_family(),
data = DryMassExp.CYBA)
Anova(betaM.1.CYBA2)
## Analysis of Deviance Table (Type II tests)
##
## Response: Moisture2
## Df Chisq Pr(>Chisq)
## Aboveground 1 7.1463 0.007512 **
## Belowground 1 2.8951 0.088847 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(betaM.1.CYBA2_glmmTMB)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Moisture2
## Chisq Df Pr(>Chisq)
## Aboveground 7.1012 1 0.007703 **
## Belowground 2.9113 1 0.087960 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(betaM.1.CYBA2_glmmTMB)
## Family: beta ( logit )
## Formula: Moisture2 ~ Aboveground + Belowground + (1 | Plot)
## Data: DryMassExp.CYBA
##
## AIC BIC logLik deviance df.resid
## -25.5 -20.7 17.7 -35.5 14
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 3.337e-10 1.827e-05
## Number of obs: 19, groups: Plot, 5
##
## Dispersion parameter for beta family (): 22.4
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.0469 0.1666 -6.285 3.28e-10 ***
## AbovegroundExposed 0.5665 0.2126 2.665 0.0077 **
## BelowgroundProtected 0.3446 0.2020 1.706 0.0880 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Answer:
Yes, there is a significant difference between
Aboveground Treatments on Moisture Content for
Cymbopetalum baillonii (X2 = 7.146, p =
0.00751). Version with RE: (X2 = 7.101, p =
0.0077).
Yes, there is a marginally significant difference between
Belowground Treatments on Moisture Content for
Cymbopetalum baillonii (X2 = 2.895, p =
0.08885). Version with RE: (X2 = 2.911, p =
0.0880).
betaM.1V2.CYBA2_glmmTMB <- glmmTMB(Moisture_NoRoots2 ~ Treatment + (1 | Plot),
family = beta_family(),
data = DryMassExp.CYBA)
summary(betaM.1V2.CYBA2_glmmTMB)
## Family: beta ( logit )
## Formula: Moisture_NoRoots2 ~ Treatment + (1 | Plot)
## Data: DryMassExp.CYBA
##
## AIC BIC logLik deviance df.resid
## -27.2 -21.5 19.6 -39.2 13
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 1.009e-10 1.004e-05
## Number of obs: 19, groups: Plot, 5
##
## Dispersion parameter for beta family (): 26
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.63956 0.23263 -2.749 0.00597 **
## TreatmentBelowground Competition 0.30796 0.32343 0.952 0.34101
## TreatmentAboveground Competition -0.04544 0.27842 -0.163 0.87036
## TreatmentFull Competition -0.67537 0.29950 -2.255 0.02414 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
DryMassExp.CYBA$Treatment <- relevel(DryMassExp.CYBA$Treatment, ref = "Belowground Competition")
betaM.1V2.CYBA2_glmmTMB <- glmmTMB(Moisture_NoRoots2 ~ Treatment + (1 | Plot),
family = beta_family(),
data = DryMassExp.CYBA)
summary(betaM.1V2.CYBA2_glmmTMB)
## Family: beta ( logit )
## Formula: Moisture_NoRoots2 ~ Treatment + (1 | Plot)
## Data: DryMassExp.CYBA
##
## AIC BIC logLik deviance df.resid
## -27.2 -21.5 19.6 -39.2 13
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 1.187e-10 1.09e-05
## Number of obs: 19, groups: Plot, 5
##
## Dispersion parameter for beta family (): 26
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.3316 0.2249 -1.475 0.140286
## TreatmentNo Competition -0.3080 0.3234 -0.952 0.341008
## TreatmentAboveground Competition -0.3534 0.2721 -1.299 0.194014
## TreatmentFull Competition -0.9833 0.2938 -3.347 0.000817 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
DryMassExp.CYBA$Treatment <- relevel(DryMassExp.CYBA$Treatment, ref = "Aboveground Competition")
betaM.1V2.CYBA2_glmmTMB <- glmmTMB(Moisture_NoRoots2 ~ Treatment + (1 | Plot),
family = beta_family(),
data = DryMassExp.CYBA)
summary(betaM.1V2.CYBA2_glmmTMB)
## Family: beta ( logit )
## Formula: Moisture_NoRoots2 ~ Treatment + (1 | Plot)
## Data: DryMassExp.CYBA
##
## AIC BIC logLik deviance df.resid
## -27.2 -21.5 19.6 -39.2 13
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 1.107e-10 1.052e-05
## Number of obs: 19, groups: Plot, 5
##
## Dispersion parameter for beta family (): 26
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.68500 0.15348 -4.463 8.08e-06 ***
## TreatmentBelowground Competition 0.35340 0.27210 1.299 0.19401
## TreatmentNo Competition 0.04544 0.27842 0.163 0.87036
## TreatmentFull Competition -0.62993 0.24313 -2.591 0.00957 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
DryMassExp.CYBA$Treatment <- relevel(DryMassExp.CYBA$Treatment, ref = "Full Competition")
betaM.1V2.CYBA2_glmmTMB <- glmmTMB(Moisture_NoRoots2 ~ Treatment + (1 | Plot),
family = beta_family(),
data = DryMassExp.CYBA)
summary(betaM.1V2.CYBA2_glmmTMB)
## Family: beta ( logit )
## Formula: Moisture_NoRoots2 ~ Treatment + (1 | Plot)
## Data: DryMassExp.CYBA
##
## AIC BIC logLik deviance df.resid
## -27.2 -21.5 19.6 -39.2 13
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 1.223e-10 1.106e-05
## Number of obs: 19, groups: Plot, 5
##
## Dispersion parameter for beta family (): 26
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.3149 0.1896 -6.936 4.04e-12 ***
## TreatmentAboveground Competition 0.6299 0.2431 2.591 0.009572 **
## TreatmentBelowground Competition 0.9833 0.2938 3.347 0.000817 ***
## TreatmentNo Competition 0.6754 0.2995 2.255 0.024136 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
betaM.1V3.CYBA2_glmmTMB <- glmmTMB(Moisture_Leaves2 ~ Treatment + (1 | Plot),
family = beta_family(),
data = DryMassExp.CYBA)
summary(betaM.1V3.CYBA2_glmmTMB)
## Family: beta ( logit )
## Formula: Moisture_Leaves2 ~ Treatment + (1 | Plot)
## Data: DryMassExp.CYBA
##
## AIC BIC logLik deviance df.resid
## -19.5 -13.8 15.8 -31.5 13
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 1.214e-09 3.485e-05
## Number of obs: 19, groups: Plot, 5
##
## Dispersion parameter for beta family (): 15.6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.4669 0.2490 -5.890 3.86e-09 ***
## TreatmentAboveground Competition 0.7115 0.3156 2.254 0.02418 *
## TreatmentBelowground Competition 1.1045 0.3793 2.912 0.00359 **
## TreatmentNo Competition 0.6562 0.3908 1.679 0.09309 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
DryMassExp.CYBA$Treatment <- relevel(DryMassExp.CYBA$Treatment, ref = "Belowground Competition")
betaM.1V3.CYBA2_glmmTMB <- glmmTMB(Moisture_Leaves2 ~ Treatment + (1 | Plot),
family = beta_family(),
data = DryMassExp.CYBA)
summary(betaM.1V3.CYBA2_glmmTMB)
## Family: beta ( logit )
## Formula: Moisture_Leaves2 ~ Treatment + (1 | Plot)
## Data: DryMassExp.CYBA
##
## AIC BIC logLik deviance df.resid
## -19.5 -13.8 15.8 -31.5 13
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 1.129e-09 3.36e-05
## Number of obs: 19, groups: Plot, 5
##
## Dispersion parameter for beta family (): 15.6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.3624 0.2871 -1.262 0.20694
## TreatmentFull Competition -1.1045 0.3793 -2.912 0.00359 **
## TreatmentAboveground Competition -0.3930 0.3481 -1.129 0.25885
## TreatmentNo Competition -0.4483 0.4176 -1.073 0.28304
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
DryMassExp.CYBA$Treatment <- relevel(DryMassExp.CYBA$Treatment, ref = "Aboveground Competition")
betaM.1V3.CYBA2_glmmTMB <- glmmTMB(Moisture_Leaves2 ~ Treatment + (1 | Plot),
family = beta_family(),
data = DryMassExp.CYBA)
summary(betaM.1V3.CYBA2_glmmTMB)
## Family: beta ( logit )
## Formula: Moisture_Leaves2 ~ Treatment + (1 | Plot)
## Data: DryMassExp.CYBA
##
## AIC BIC logLik deviance df.resid
## -19.5 -13.8 15.8 -31.5 13
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 8.95e-10 2.992e-05
## Number of obs: 19, groups: Plot, 5
##
## Dispersion parameter for beta family (): 15.6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.75540 0.19741 -3.826 0.00013 ***
## TreatmentBelowground Competition 0.39303 0.34808 1.129 0.25885
## TreatmentFull Competition -0.71147 0.31561 -2.254 0.02418 *
## TreatmentNo Competition -0.05523 0.36129 -0.153 0.87850
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
DryMassExp.CYBA$Treatment <- relevel(DryMassExp.CYBA$Treatment, ref = "No Competition")
betaM.1V3.CYBA2_glmmTMB <- glmmTMB(Moisture_Leaves2 ~ Treatment + (1 | Plot),
family = beta_family(),
data = DryMassExp.CYBA)
summary(betaM.1V3.CYBA2_glmmTMB)
## Family: beta ( logit )
## Formula: Moisture_Leaves2 ~ Treatment + (1 | Plot)
## Data: DryMassExp.CYBA
##
## AIC BIC logLik deviance df.resid
## -19.5 -13.8 15.8 -31.5 13
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot (Intercept) 8.656e-10 2.942e-05
## Number of obs: 19, groups: Plot, 5
##
## Dispersion parameter for beta family (): 15.6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.81063 0.30362 -2.670 0.00759 **
## TreatmentAboveground Competition 0.05523 0.36129 0.153 0.87850
## TreatmentBelowground Competition 0.44826 0.41756 1.073 0.28304
## TreatmentFull Competition -0.65624 0.39077 -1.679 0.09309 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Answer:
T7B$Total <- 10
T7B$Establishment <- T7B$Count/ 10 * 100
tgc <- summarySE(T7B, measurevar="Establishment",
groupvars=c("Species"), na.rm = TRUE)
levels(tgc$Species) <- list("SASA" = "Sapindus saponaria",
"CYBA" = "Cymbopetalum baillonii")
GG_Est_Spp <- ggplot(tgc, aes(x=Species, y=Establishment, fill=Species)) +
geom_bar(stat = "identity", color = "black", alpha = .5) +
geom_errorbar(aes(ymin=Establishment-se, ymax=Establishment+se), width=.1) +
geom_text(label = c("a", "b"), aes(y = Establishment+se, x = Species),vjust = -0.5, size = 5, family = "serif") +
scale_fill_manual(values=c("darkslategray3", "coral2")) +
ylab("Seedling \nEstablishment (%)") +
ylim(0, 42) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
text = element_text(family = "serif"),
axis.title = element_text(size = rel(1)),
axis.text = element_text(size = rel(0.9)),
axis.ticks = element_line(color = "black", size = 0.8),
panel.border = element_rect(color = "black", fill = NA, size = 0.8),
strip.background = element_rect(color = "black", fill = "white", size = 0.8),
panel.grid = element_blank())
GG_Est_Spp
png("~/Desktop/Universidad Javeriana/Manuscripts/Fern Interference/Biotropica/Prototype Figures/Establishment x Species.png", width = 70, height = 80, units = 'mm', res = 300)
GG_Est_Spp
dev.off()
## quartz_off_screen
## 2
tgc <- summarySE(T7B, measurevar="Establishment",
groupvars=c("Aboveground" , "Species"), na.rm = TRUE)
tgc$Species_Aboveground <- paste(tgc$Species, tgc$Aboveground, sep = "_")
levels(tgc$Aboveground) <- list("AG C." = "Covered",
"No AG C." = "Exposed")
tgc$Species <- relevel(tgc$Species, ref = "Sapindus saponaria")
levels(tgc$Species) <- list("Sapindus" = "Sapindus saponaria",
"Cymbopetalum" = "Cymbopetalum baillonii")
tgc$Aboveground <- relevel(tgc$Aboveground, ref = "AG C.")
GG_Est_AG <- ggplot(tgc, aes(x=Aboveground, y=Establishment, fill=Species_Aboveground)) +
geom_bar(stat = "identity", color = "black", alpha = .5) +
geom_errorbar(aes(ymin=Establishment-se, ymax=Establishment+se), width=.1) +
geom_text(label = c("a", "a", "b", "a"),
aes(y = Establishment+se, x = Aboveground),vjust = -0.5, size = 5, family = "serif") +
scale_fill_manual(values = list("paleturquoise", "paleturquoise4",
"lightsalmon", "lightsalmon4")) +
ylab("Seedling \nEstablishment (%)") +
xlab("Aboveground Treatment") +
ylim(0, 55) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
text = element_text(family = "serif"),
axis.title = element_text(size = rel(1)),
axis.text = element_text(size = rel(0.9)),
axis.ticks = element_line(color = "black", size = 0.8),
panel.border = element_rect(color = "black", fill = NA, size = 0.8),
strip.background = element_rect(color = "black", fill = "white", size = 0.8),
panel.grid = element_blank()) +
facet_wrap(~Species, scales = "free", strip.position = "top")
GG_Est_AG
png("~/Desktop/Universidad Javeriana/Manuscripts/Fern Interference/Biotropica/Prototype Figures/Establishment x AG.png", width = 70, height = 80, units = 'mm', res = 300)
GG_Est_AG
dev.off()
## quartz_off_screen
## 2
tgc <- summarySE(T7B, measurevar="Establishment",
groupvars=c("Belowground", "Species"), na.rm = TRUE)
tgc$Species_Belowground <- paste(tgc$Species, tgc$Belowground, sep = "_")
levels(tgc$Belowground) <- list("BG C." = "Unprotected",
"No BG C." = "Protected")
tgc$Species <- relevel(tgc$Species, ref = "Sapindus saponaria")
levels(tgc$Species) <- list("Sapindus" = "Sapindus saponaria",
"Cymbopetalum" = "Cymbopetalum baillonii")
tgc$Belowground <- relevel(tgc$Belowground, ref = "BG C.")
GG_Est_BG <- ggplot(tgc, aes(x=Belowground, y=Establishment, fill=Species_Belowground)) +
geom_bar(stat = "identity", color = "black", alpha = .5) +
geom_errorbar(aes(ymin=Establishment-se, ymax=Establishment+se), width=.1) +
geom_text(label = c("a", "a", "a", "a"), aes(y = Establishment+se, x = Belowground),vjust = -0.5, size = 5, family = "serif") +
scale_fill_manual(values = list("skyblue4", "skyblue",
"mistyrose4", "mistyrose")) +
ylab("Seedling \nEstablishment (%)") +
xlab("Belowground Treatment") +
ylim(0, 50) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
text = element_text(family = "serif"),
axis.title = element_text(size = rel(1)),
axis.text = element_text(size = rel(0.9)),
axis.ticks = element_line(color = "black", size = 0.8),
panel.border = element_rect(color = "black", fill = NA, size = 0.8),
strip.background = element_rect(color = "black", fill = "white", size = 0.8),
panel.grid = element_blank()) +
facet_wrap(~Species, scales = "free", strip.position = "top")
GG_Est_BG
png("~/Desktop/Universidad Javeriana/Manuscripts/Fern Interference/Biotropica/Prototype Figures/Establishment x BG.png", width = 70, height = 80, units = 'mm', res = 300)
GG_Est_BG
dev.off()
## quartz_off_screen
## 2
tgc <- summarySE(FernExp.G2, measurevar="MaxLeaves",
groupvars=c("Species"), na.rm = TRUE)
levels(tgc$Species) <- list("SASA" = "Sapindus saponaria",
"CYBA" = "Cymbopetalum baillonii")
GG_LeafCount_Spp <- ggplot(tgc, aes(x=Species, y=MaxLeaves, fill=Species)) +
geom_bar(stat = "identity", color = "black", alpha = .5) +
geom_errorbar(aes(ymin=MaxLeaves-se, ymax=MaxLeaves+se), width=.1) +
geom_text(label = c("b", "a"), aes(y = MaxLeaves+se, x = Species),vjust = -0.5, size = 5, family = "serif") +
scale_fill_manual(values=c("darkslategray3", "coral2")) +
ylab("\nHighest Leaf Count") +
ylim(0, 10) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
text = element_text(family = "serif"),
axis.title = element_text(size = rel(1)),
axis.text = element_text(size = rel(0.9)), axis.ticks = element_line(color = "black", size = 0.8),
panel.border = element_rect(color = "black", fill = NA, size = 0.8),
strip.background = element_rect(color = "black", fill = "white", size = 0.8),
panel.grid = element_blank())
GG_LeafCount_Spp
png("~/Desktop/Universidad Javeriana/Manuscripts/Fern Interference/Biotropica/Prototype Figures/Leaf Count x Species.png", width = 70, height = 80, units = 'mm', res = 300)
GG_LeafCount_Spp
dev.off()
## quartz_off_screen
## 2
tgc
## Species N MaxLeaves sd se ci
## 1 CYBA 78 4.423077 5.616578 0.6359523 1.266343
## 2 SASA 47 8.085106 7.606718 1.1095539 2.233416
(4.42308 - 8.08511) / 4.42308 * 100
## [1] -82.79366
tgc <- summarySE(FernExp.G2, measurevar="MaxLeaves",
groupvars=c("Aboveground", "Species"), na.rm = TRUE)
tgc$Species_Aboveground <- paste(tgc$Species, tgc$Aboveground, sep = "_")
levels(tgc$Aboveground) <- list("AG C." = "Covered",
"No AG C." = "Exposed")
tgc$Species <- relevel(tgc$Species, ref = "Sapindus saponaria")
levels(tgc$Species) <- list("Sapindus" = "Sapindus saponaria",
"Cymbopetalum" = "Cymbopetalum baillonii")
tgc$Aboveground <- relevel(tgc$Aboveground, ref = "AG C.")
GG_LeafCount_AG <- ggplot(tgc, aes(x=Aboveground, y=MaxLeaves, fill=Species_Aboveground)) +
geom_bar(stat = "identity", color = "black", alpha = .5) +
geom_errorbar(aes(ymin=MaxLeaves-se, ymax=MaxLeaves+se), width=.1) +
geom_text(label = c("a", "b", "a", "a"), aes(y = MaxLeaves+se, x = Aboveground),vjust = -0.5, size = 5, family = "serif") +
scale_fill_manual(values = list("paleturquoise", "paleturquoise4",
"lightsalmon", "lightsalmon4")) +
ylab("\nHighest Leaf Count") +
ylim(0, 15) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
text = element_text(family = "serif"),
axis.title = element_text(size = rel(1)),
axis.text = element_text(size = rel(0.9)), axis.ticks = element_line(color = "black", size = 0.8),
panel.border = element_rect(color = "black", fill = NA, size = 0.8),
strip.background = element_rect(color = "black", fill = "white", size = 0.8),
panel.grid = element_blank()) +
facet_wrap(~Species, scales = "free", strip.position = "top")
GG_LeafCount_AG
png("~/Desktop/Universidad Javeriana/Manuscripts/Fern Interference/Biotropica/Prototype Figures/Leaf Count x AG.png", width = 70, height = 80, units = 'mm', res = 300)
GG_LeafCount_AG
dev.off()
## quartz_off_screen
## 2
tgc <- summarySE(FernExp.G2, measurevar="MaxLeaves",
groupvars=c("Belowground", "Species"), na.rm = TRUE)
tgc$Species_Belowground <- paste(tgc$Species, tgc$Belowground, sep = "_")
levels(tgc$Belowground) <- list("BG C." = "Unprotected",
"No BG C." = "Protected")
tgc$Species <- relevel(tgc$Species, ref = "Sapindus saponaria")
levels(tgc$Species) <- list("Sapindus" = "Sapindus saponaria",
"Cymbopetalum" = "Cymbopetalum baillonii")
tgc$Belowground <- relevel(tgc$Belowground, ref = "BG C.")
GG_LeafCount_BG <- ggplot(tgc, aes(x=Belowground, y=MaxLeaves, fill=Species_Belowground)) +
geom_bar(stat = "identity", color = "black", alpha = .5) +
geom_errorbar(aes(ymin=MaxLeaves-se, ymax=MaxLeaves+se), width=.1) +
geom_text(label = c("a", "a", "a", "a"), aes(y = MaxLeaves+se, x = Belowground),vjust = -0.5, size = 5, family = "serif") +
scale_fill_manual(values = list("skyblue4", "skyblue","mistyrose4", "mistyrose")) +
ylab("\nHighest Leaf Count") +
ylim(0, 10.5) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
text = element_text(family = "serif"),
axis.title = element_text(size = rel(1)),
axis.text = element_text(size = rel(0.9)), axis.ticks = element_line(color = "black", size = 0.8),
panel.border = element_rect(color = "black", fill = NA, size = 0.8),
strip.background = element_rect(color = "black", fill = "white", size = 0.8),
panel.grid = element_blank()) +
facet_wrap(~Species, scales = "free", strip.position = "top")
GG_LeafCount_BG
png("~/Desktop/Universidad Javeriana/Manuscripts/Fern Interference/Biotropica/Prototype Figures/Leaf Count x BG.png", width = 70, height = 80, units = 'mm', res = 300)
GG_LeafCount_BG
dev.off()
## quartz_off_screen
## 2
tgc <- summarySE(WetMassExp, measurevar="Height..cm.",
groupvars=c("Species"), na.rm = TRUE)
levels(tgc$Species) <- list("SASA" = "Sapindus saponaria",
"CYBA" = "Cymbopetalum baillonii")
GG_Height_Spp <- ggplot(tgc, aes(x=Species, y=Height..cm., fill=Species)) +
geom_bar(stat = "identity", color = "black", alpha = .5) +
geom_errorbar(aes(ymin=Height..cm.-se, ymax=Height..cm.+se), width=.1) +
geom_text(label = c("a", "a"), aes(y = Height..cm.+se, x = Species),vjust = -0.5, size = 5, family = "serif") +
scale_fill_manual(values=c("darkslategray3", "coral2")) +
ylab("Seedling \nHeight (cm)") +
ylim(0, 30) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
text = element_text(family = "serif"),
axis.title = element_text(size = rel(1)),
axis.text = element_text(size = rel(0.9)), axis.ticks = element_line(color = "black", size = 0.8),
panel.border = element_rect(color = "black", fill = NA, size = 0.8),
strip.background = element_rect(color = "black", fill = "white", size = 0.8),
panel.grid = element_blank())
GG_Height_Spp
png("~/Desktop/Universidad Javeriana/Manuscripts/Fern Interference/Biotropica/Prototype Figures/Height x Species.png", width = 70, height = 80, units = 'mm', res = 300)
GG_Height_Spp
dev.off()
## quartz_off_screen
## 2
tgc <- summarySE(WetMassExp, measurevar="Height..cm.",
groupvars=c("Aboveground", "Species"), na.rm = TRUE)
tgc$Species_Aboveground <- paste(tgc$Species, tgc$Aboveground, sep = "_")
levels(tgc$Aboveground) <- list("AG C." = "Covered",
"No AG C." = "Exposed")
tgc$Species <- relevel(tgc$Species, ref = "Sapindus saponaria")
levels(tgc$Species) <- list("Sapindus" = "Sapindus saponaria",
"Cymbopetalum" = "Cymbopetalum baillonii")
tgc$Aboveground <- relevel(tgc$Aboveground, ref = "AG C.")
GG_Height_AG <- ggplot(tgc, aes(x=Aboveground, y=Height..cm., fill=Species_Aboveground)) +
geom_bar(stat = "identity", color = "black", alpha = .5) +
geom_errorbar(aes(ymin=Height..cm.-se, ymax=Height..cm.+se), width=.1) +
geom_text(label = c("a", "b","a", "a"), aes(y = Height..cm.+se, x = Aboveground),vjust = -0.5, size = 5, family = "serif") +
scale_fill_manual(values = list("paleturquoise", "paleturquoise4","lightsalmon", "lightsalmon4")) +
ylab("Seedling \nHeight (cm)") +
ylim(0, 43) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
text = element_text(family = "serif"),
axis.title = element_text(size = rel(1)),
axis.text = element_text(size = rel(0.9)), axis.ticks = element_line(color = "black", size = 0.8),
panel.border = element_rect(color = "black", fill = NA, size = 0.8),
strip.background = element_rect(color = "black", fill = "white", size = 0.8),
panel.grid = element_blank()) +
facet_wrap(~Species, scales = "free", strip.position = "top")
GG_Height_AG
png("~/Desktop/Universidad Javeriana/Manuscripts/Fern Interference/Biotropica/Prototype Figures/Height x AG.png", width = 70, height = 80, units = 'mm', res = 300)
GG_Height_AG
dev.off()
## quartz_off_screen
## 2
tgc
## Aboveground Species N Height..cm. sd se ci
## 1 AG C. Cymbopetalum 13 19.00 5.730038 1.589227 3.462627
## 2 AG C. Sapindus 15 14.40 6.009516 1.551650 3.327959
## 3 No AG C. Cymbopetalum 6 20.50 10.502381 4.287579 11.021573
## 4 No AG C. Sapindus 16 32.25 30.708305 7.677076 16.363301
## Species_Aboveground
## 1 Cymbopetalum baillonii_Covered
## 2 Sapindus saponaria_Covered
## 3 Cymbopetalum baillonii_Exposed
## 4 Sapindus saponaria_Exposed
(32.25 / 14.40) / 32.25 * 100
## [1] 6.944444
tgc <- summarySE(WetMassExp, measurevar="Height..cm.",
groupvars=c("Belowground", "Species"), na.rm = TRUE)
tgc$Species_Belowground <- paste(tgc$Species, tgc$Belowground, sep = "_")
levels(tgc$Belowground) <- list("BG C." = "Unprotected",
"No BG C." = "Protected")
tgc$Species <- relevel(tgc$Species, ref = "Sapindus saponaria")
levels(tgc$Species) <- list("Sapindus" = "Sapindus saponaria",
"Cymbopetalum" = "Cymbopetalum baillonii")
tgc$Belowground <- relevel(tgc$Belowground, ref = "BG C.")
GG_Height_BG <- ggplot(tgc, aes(x=Belowground, y=Height..cm., fill=Species_Belowground)) +
geom_bar(stat = "identity", color = "black", alpha = .5) +
geom_errorbar(aes(ymin=Height..cm.-se, ymax=Height..cm.+se), width=.1) +
geom_text(label = c("a", "a", "a", "a"), aes(y = Height..cm.+se, x = Belowground),vjust = -0.5, size = 5, family = "serif") +
scale_fill_manual(values = list("skyblue4", "skyblue","mistyrose4", "mistyrose")) +
ylab("Seedling \nHeight (cm)") +
ylim(0, 36) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
text = element_text(family = "serif"),
axis.title = element_text(size = rel(1)),
axis.text = element_text(size = rel(0.9)), axis.ticks = element_line(color = "black", size = 0.8),
panel.border = element_rect(color = "black", fill = NA, size = 0.8),
strip.background = element_rect(color = "black", fill = "white", size = 0.8),
panel.grid = element_blank()) +
facet_wrap(~Species, scales = "free", strip.position = "top")
GG_Height_BG
png("~/Desktop/Universidad Javeriana/Manuscripts/Fern Interference/Biotropica/Prototype Figures/Height x BG.png", width = 70, height = 80, units = 'mm', res = 300)
GG_Height_BG
dev.off()
## quartz_off_screen
## 2
tgc <- summarySE(WetMassExp, measurevar="Diameter.Base..cm.",
groupvars=c("Species"), na.rm = TRUE)
levels(tgc$Species) <- list("SASA" = "Sapindus saponaria",
"CYBA" = "Cymbopetalum baillonii")
GG_Diam_Spp <- ggplot(tgc, aes(x=Species, y=Diameter.Base..cm., fill=Species)) +
geom_bar(stat = "identity", color = "black", alpha = .5) +
geom_errorbar(aes(ymin=Diameter.Base..cm.-se, ymax=Diameter.Base..cm.+se), width=.1) +
geom_text(label = c("a", "a"), aes(y = Diameter.Base..cm.+se, x = Species),vjust = -0.5, size = 5, family = "serif") +
scale_fill_manual(values=c("darkslategray3", "coral2")) +
ylab("Base \nDiameter (cm)") +
ylim(0, 1.25) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
text = element_text(family = "serif"),
axis.title = element_text(size = rel(1)),
axis.text = element_text(size = rel(0.9)), axis.ticks = element_line(color = "black", size = 0.8),
panel.border = element_rect(color = "black", fill = NA, size = 0.8),
strip.background = element_rect(color = "black", fill = "white", size = 0.8),
panel.grid = element_blank())
GG_Diam_Spp
png("~/Desktop/Universidad Javeriana/Manuscripts/Fern Interference/Biotropica/Prototype Figures/Diameter x Species.png", width = 70, height = 80, units = 'mm', res = 300)
GG_Diam_Spp
dev.off()
## quartz_off_screen
## 2
tgc <- summarySE(WetMassExp, measurevar="Diameter.Base..cm.",
groupvars=c("Aboveground", "Species"), na.rm = TRUE)
tgc$Species_Aboveground <- paste(tgc$Species, tgc$Aboveground, sep = "_")
levels(tgc$Aboveground) <- list("AG C." = "Covered",
"No AG C." = "Exposed")
tgc$Species <- relevel(tgc$Species, ref = "Sapindus saponaria")
levels(tgc$Species) <- list("Sapindus" = "Sapindus saponaria",
"Cymbopetalum" = "Cymbopetalum baillonii")
tgc$Aboveground <- relevel(tgc$Aboveground, ref = "AG C.")
GG_Diam_AG <- ggplot(tgc, aes(x=Aboveground, y=Diameter.Base..cm., fill=Species_Aboveground)) +
geom_bar(stat = "identity", color = "black", alpha = .5) +
geom_errorbar(aes(ymin=Diameter.Base..cm.-se, ymax=Diameter.Base..cm.+se), width=.1) +
geom_text(label = c("b", "b", "a", "a"), aes(y = Diameter.Base..cm.+se, x = Aboveground),vjust = -0.5, size = 5, family = "serif") +
scale_fill_manual(values = list("paleturquoise", "paleturquoise4","lightsalmon", "lightsalmon4")) +
ylab("\nBase Diameter (cm)") +
ylim(0, 2.5) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
text = element_text(family = "serif"),
axis.title = element_text(size = rel(1)),
axis.text = element_text(size = rel(0.9)), axis.ticks = element_line(color = "black", size = 0.8),
panel.border = element_rect(color = "black", fill = NA, size = 0.8),
strip.background = element_rect(color = "black", fill = "white", size = 0.8),
panel.grid = element_blank()) +
facet_wrap(~Species, scales = "free", strip.position = "top")
GG_Diam_AG
png("~/Desktop/Universidad Javeriana/Manuscripts/Fern Interference/Biotropica/Prototype Figures/Diameter x AG.png", width = 70, height = 80, units = 'mm', res = 300)
GG_Diam_AG
dev.off()
## quartz_off_screen
## 2
tgc
## Aboveground Species N Diameter.Base..cm. sd se
## 1 AG C. Cymbopetalum 13 0.3538462 0.07762500 0.02152930
## 2 AG C. Sapindus 15 0.3066667 0.08837151 0.02281743
## 3 No AG C. Cymbopetalum 6 0.5500000 0.33911650 0.13844373
## 4 No AG C. Sapindus 16 1.3875000 2.38463834 0.59615958
## ci Species_Aboveground
## 1 0.04690832 Cymbopetalum baillonii_Covered
## 2 0.04893851 Sapindus saponaria_Covered
## 3 0.35588094 Cymbopetalum baillonii_Exposed
## 4 1.27068408 Sapindus saponaria_Exposed
(0.306667 - 1.387500) / 0.306667 * 100
## [1] -352.4452
(0.353846 - 0.550000) / 0.306667 * 100
## [1] -63.96319
tgc <- summarySE(WetMassExp, measurevar="Diameter.Base..cm.",
groupvars=c("Belowground", "Species"), na.rm = TRUE)
tgc$Species_Belowground <- paste(tgc$Species, tgc$Belowground, sep = "_")
levels(tgc$Belowground) <- list("BG C." = "Unprotected",
"No BG C." = "Protected")
tgc$Species <- relevel(tgc$Species, ref = "Sapindus saponaria")
levels(tgc$Species) <- list("Sapindus" = "Sapindus saponaria",
"Cymbopetalum" = "Cymbopetalum baillonii")
tgc$Belowground <- relevel(tgc$Belowground, ref = "BG C.")
GG_Diam_BG <- ggplot(tgc, aes(x=Belowground, y=Diameter.Base..cm., fill=Species_Belowground)) +
geom_bar(stat = "identity", color = "black", alpha = .5) +
geom_errorbar(aes(ymin=Diameter.Base..cm.-se, ymax=Diameter.Base..cm.+se), width=.1) +
geom_text(label = c("a", " b·", "a", " a·"), aes(y = Diameter.Base..cm.+se, x = Belowground),vjust = -0.5, size = 5, family = "serif") +
scale_fill_manual(values = list("skyblue4", "skyblue","mistyrose4", "mistyrose")) +
ylab("\nBase Diameter (cm)") +
ylim(0, 2.5) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
text = element_text(family = "serif"),
axis.title = element_text(size = rel(1)),
axis.text = element_text(size = rel(0.9)), axis.ticks = element_line(color = "black", size = 0.8),
panel.border = element_rect(color = "black", fill = NA, size = 0.8),
strip.background = element_rect(color = "black", fill = "white", size = 0.8),
panel.grid = element_blank()) +
facet_wrap(~Species, scales = "free", strip.position = "top")
GG_Diam_BG
png("~/Desktop/Universidad Javeriana/Manuscripts/Fern Interference/Biotropica/Prototype Figures/Diameter x BG.png", width = 70, height = 80, units = 'mm', res = 300)
GG_Diam_BG
dev.off()
## quartz_off_screen
## 2
tgc
## Belowground Species N Diameter.Base..cm. sd se
## 1 No BG C. Cymbopetalum 10 0.3800000 0.09189366 0.02905933
## 2 No BG C. Sapindus 17 0.5176471 0.32061522 0.07776061
## 3 BG C. Cymbopetalum 9 0.4555556 0.29627315 0.09875772
## 4 BG C. Sapindus 14 1.2857143 2.60586827 0.69644759
## ci Species_Belowground
## 1 0.06573676 Cymbopetalum baillonii_Protected
## 2 0.16484513 Sapindus saponaria_Protected
## 3 0.22773570 Cymbopetalum baillonii_Unprotected
## 4 1.50458355 Sapindus saponaria_Unprotected
(0.517647 - 1.285714) / 0.517647 * 100
## [1] -148.3766
(1.285714 - 0.517647) / 1.285714 * 100
## [1] 59.73856
tgc <- summarySE(DryMassExp, measurevar="TotalDryMass",
groupvars=c("Species"), na.rm = TRUE)
levels(tgc$Species) <- list("SASA" = "Sapindus saponaria",
"CYBA" = "Cymbopetalum baillonii")
GG_DMass_Spp <- ggplot(tgc, aes(x=Species, y=TotalDryMass, fill=Species)) +
geom_bar(stat = "identity", color = "black", alpha = .5) +
geom_errorbar(aes(ymin=TotalDryMass-se, ymax=TotalDryMass+se), width=.1) +
geom_text(label = c("b·", "a·"), aes(y = TotalDryMass+se, x = Species),vjust = -0.5, size = 5, family = "serif") +
scale_fill_manual(values=c("darkslategray3", "coral2")) +
ylab("\nDry Total Mass (g)") +
ylim(0, 21) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
text = element_text(family = "serif"),
axis.title = element_text(size = rel(1)),
axis.text = element_text(size = rel(0.9)), axis.ticks = element_line(color = "black", size = 0.8),
panel.border = element_rect(color = "black", fill = NA, size = 0.8),
strip.background = element_rect(color = "black", fill = "white", size = 0.8),
panel.grid = element_blank())
GG_DMass_Spp
png("~/Desktop/Universidad Javeriana/Manuscripts/Fern Interference/Biotropica/Prototype Figures/Dry Mass x Species.png", width = 70, height = 80, units = 'mm', res = 300)
GG_DMass_Spp
dev.off()
## quartz_off_screen
## 2
tgc
## Species N TotalDryMass sd se ci
## 1 CYBA 19 1.732632 2.636445 0.604842 1.270726
## 2 SASA 31 12.000968 38.164386 6.854526 13.998810
(1.73263 - 12.00097) / 1.73263 * 100
## [1] -592.6447
tgc <- summarySE(DryMassExp, measurevar="TotalDryMass",
groupvars=c("Aboveground", "Species"), na.rm = TRUE)
tgc$Species_Aboveground <- paste(tgc$Species, tgc$Aboveground, sep = "_")
levels(tgc$Aboveground) <- list("AG C." = "Covered",
"No AG C." = "Exposed")
tgc$Species <- relevel(tgc$Species, ref = "Sapindus saponaria")
levels(tgc$Species) <- list("Sapindus" = "Sapindus saponaria",
"Cymbopetalum" = "Cymbopetalum baillonii")
tgc$Aboveground <- relevel(tgc$Aboveground, ref = "AG C.")
GG_DMass_AG <- ggplot(tgc, aes(x=Aboveground, y=TotalDryMass, fill=Species_Aboveground)) +
geom_bar(stat = "identity", color = "black", alpha = .5) +
geom_errorbar(aes(ymin=TotalDryMass-se, ymax=TotalDryMass+se), width=.1) +
geom_text(label = c("b", "b", "a", "a"), aes(y = TotalDryMass+se, x = Aboveground),vjust = -0.5, size = 5, family = "serif") +
scale_fill_manual(values = list("paleturquoise", "paleturquoise4","lightsalmon", "lightsalmon4")) +
ylab("\nDry Seedling Mass (g)") +
ylim(0, 40) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
text = element_text(family = "serif"),
axis.title = element_text(size = rel(1)),
axis.text = element_text(size = rel(0.9)), axis.ticks = element_line(color = "black", size = 0.8),
panel.border = element_rect(color = "black", fill = NA, size = 0.8),
strip.background = element_rect(color = "black", fill = "white", size = 0.8),
panel.grid = element_blank()) +
facet_wrap(~Species, scales = "free", strip.position = "top")
GG_DMass_AG
png("~/Desktop/Universidad Javeriana/Manuscripts/Fern Interference/Biotropica/Prototype Figures/Dry Mass x AG.png", width = 70, height = 80, units = 'mm', res = 300)
GG_DMass_AG
dev.off()
## quartz_off_screen
## 2
tgc
## Aboveground Species N TotalDryMass sd se ci
## 1 AG C. Cymbopetalum 13 0.8407692 0.6141588 0.1703370 0.3711324
## 2 AG C. Sapindus 15 0.7586667 0.6247384 0.1613068 0.3459686
## 3 No AG C. Cymbopetalum 6 3.6650000 4.1915188 1.7111804 4.3987292
## 4 No AG C. Sapindus 16 22.5406250 51.6507227 12.9126807 27.5227273
## Species_Aboveground
## 1 Cymbopetalum baillonii_Covered
## 2 Sapindus saponaria_Covered
## 3 Cymbopetalum baillonii_Exposed
## 4 Sapindus saponaria_Exposed
(0.758667 - 22.540625) / 0.758667 * 100
## [1] -2871.083
(0.840769 - 3.665000) / 0.840769 * 100
## [1] -335.9105
tgc <- summarySE(DryMassExp, measurevar="TotalDryMass",
groupvars=c("Belowground", "Species"), na.rm = TRUE)
tgc$Species_Belowground <- paste(tgc$Species, tgc$Belowground, sep = "_")
levels(tgc$Belowground) <- list("BG C." = "Unprotected",
"No BG C." = "Protected")
tgc$Species <- relevel(tgc$Species, ref = "Sapindus saponaria")
levels(tgc$Species) <- list("Sapindus" = "Sapindus saponaria",
"Cymbopetalum" = "Cymbopetalum baillonii")
tgc$Belowground <- relevel(tgc$Belowground, ref = "BG C.")
GG_DMass_BG <- ggplot(tgc, aes(x=Belowground, y=TotalDryMass, fill=Species_Belowground)) +
geom_bar(stat = "identity", color = "black", alpha = .5) +
geom_errorbar(aes(ymin=TotalDryMass-se, ymax=TotalDryMass+se), width=.1) +
geom_text(label = c("a", "a", "a", "a"), aes(y = TotalDryMass+se, x = Belowground),vjust = -0.5, size = 5, family = "serif") +
scale_fill_manual(values = list("skyblue4", "skyblue","mistyrose4", "mistyrose")) +
ylab("\nDry Seedling Mass (g)") +
ylim(0, 36) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
text = element_text(family = "serif"),
axis.title = element_text(size = rel(1)),
axis.text = element_text(size = rel(0.9)), axis.ticks = element_line(color = "black", size = 0.8),
panel.border = element_rect(color = "black", fill = NA, size = 0.8),
strip.background = element_rect(color = "black", fill = "white", size = 0.8),
panel.grid = element_blank()) +
facet_wrap(~Species, scales = "free", strip.position = "top")
GG_DMass_BG
png("~/Desktop/Universidad Javeriana/Manuscripts/Fern Interference/Biotropica/Prototype Figures/Dry Mass x BG.png", width = 70, height = 80, units = 'mm', res = 300)
GG_DMass_BG
dev.off()
## quartz_off_screen
## 2
tgc <- summarySE(DryMassExp, measurevar="Ratio",
groupvars=c("Species"), na.rm = TRUE)
levels(tgc$Species) <- list("SASA" = "Sapindus saponaria",
"CYBA" = "Cymbopetalum baillonii")
GG_Drati_Spp <- ggplot(tgc, aes(x=Species, y=Ratio, fill=Species)) +
geom_bar(stat = "identity", color = "black", alpha = .5) +
geom_errorbar(aes(ymin=Ratio-se, ymax=Ratio+se), width=.1) +
geom_text(label = c("b", "a"), aes(y = Ratio+se, x = Species),vjust = -0.5, size = 5, family = "serif") +
scale_fill_manual(values=c("darkslategray3", "coral2")) +
ylab("\nDry R/S Ratio") +
ylim(0, 1.3) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
text = element_text(family = "serif"),
axis.title = element_text(size = rel(1)),
axis.text = element_text(size = rel(0.9)), axis.ticks = element_line(color = "black", size = 0.8),
panel.border = element_rect(color = "black", fill = NA, size = 0.8),
strip.background = element_rect(color = "black", fill = "white", size = 0.8),
panel.grid = element_blank())
GG_Drati_Spp
png("~/Desktop/Universidad Javeriana/Manuscripts/Fern Interference/Biotropica/Prototype Figures/Dry Ratio x Species.png", width = 70, height = 80, units = 'mm', res = 300)
GG_Drati_Spp
dev.off()
## quartz_off_screen
## 2
tgc
## Species N Ratio sd se ci
## 1 CYBA 19 0.6910584 0.3721360 0.08537385 0.1793638
## 2 SASA 31 1.1237638 0.4740003 0.08513296 0.1738647
(0.691058 - 1.123764) / 0.691058 *100
## [1] -62.615
tgc <- summarySE(DryMassExp, measurevar="Ratio",
groupvars=c("Aboveground", "Species"), na.rm = TRUE)
tgc$Species_Aboveground <- paste(tgc$Species, tgc$Aboveground, sep = "_")
levels(tgc$Aboveground) <- list("AG C." = "Covered",
"No AG C." = "Exposed")
tgc$Species <- relevel(tgc$Species, ref = "Sapindus saponaria")
levels(tgc$Species) <- list("Sapindus" = "Sapindus saponaria",
"Cymbopetalum" = "Cymbopetalum baillonii")
tgc$Aboveground <- relevel(tgc$Aboveground, ref = "AG C.")
GG_Drati_AG <- ggplot(tgc, aes(x=Aboveground, y=Ratio, fill=Species_Aboveground)) +
geom_bar(stat = "identity", color = "black", alpha = .5) +
geom_errorbar(aes(ymin=Ratio-se, ymax=Ratio+se), width=.1) +
geom_text(label = c("a", "a", "a", "a"), aes(y = Ratio+se, x = Aboveground),vjust = -0.5, size = 5, family = "serif") +
scale_fill_manual(values = list("paleturquoise", "paleturquoise4","lightsalmon", "lightsalmon4")) +
ylab("\nDry R/S Ratio") +
ylim(0, 1.5) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
text = element_text(family = "serif"),
axis.title = element_text(size = rel(1)),
axis.text = element_text(size = rel(0.9)), axis.ticks = element_line(color = "black", size = 0.8),
panel.border = element_rect(color = "black", fill = NA, size = 0.8),
strip.background = element_rect(color = "black", fill = "white", size = 0.8),
panel.grid = element_blank()) +
facet_wrap(~Species, scales = "free", strip.position = "top")
GG_Drati_AG
png("~/Desktop/Universidad Javeriana/Manuscripts/Fern Interference/Biotropica/Prototype Figures/Dry Ratio x AG.png", width = 70, height = 80, units = 'mm', res = 300)
GG_Drati_AG
dev.off()
## quartz_off_screen
## 2
tgc <- summarySE(DryMassExp, measurevar="Ratio",
groupvars=c("Belowground", "Species"), na.rm = TRUE)
tgc$Species_Belowground <- paste(tgc$Species, tgc$Belowground, sep = "_")
levels(tgc$Belowground) <- list("BG C." = "Unprotected",
"No BG C." = "Protected")
tgc$Species <- relevel(tgc$Species, ref = "Sapindus saponaria")
levels(tgc$Species) <- list("Sapindus" = "Sapindus saponaria",
"Cymbopetalum" = "Cymbopetalum baillonii")
tgc$Belowground <- relevel(tgc$Belowground, ref = "BG C.")
GG_Drati_BG <- ggplot(tgc, aes(x=Belowground, y=Ratio, fill=Species_Belowground)) +
geom_bar(stat = "identity", color = "black", alpha = .5) +
geom_errorbar(aes(ymin=Ratio-se, ymax=Ratio+se), width=.1) +
geom_text(label = c("a", "a", "a", "a"), aes(y = Ratio+se, x = Belowground),vjust = -0.5, size = 5, family = "serif") +
scale_fill_manual(values = list("skyblue4", "skyblue","mistyrose4", "mistyrose")) +
ylab("\nDry R/S Ratio") +
ylim(0, 1.5) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
text = element_text(family = "serif"),
axis.title = element_text(size = rel(1)),
axis.text = element_text(size = rel(0.9)), axis.ticks = element_line(color = "black", size = 0.8),
panel.border = element_rect(color = "black", fill = NA, size = 0.8),
strip.background = element_rect(color = "black", fill = "white", size = 0.8),
panel.grid = element_blank()) +
facet_wrap(~Species, scales = "free", strip.position = "top")
GG_Drati_BG
png("~/Desktop/Universidad Javeriana/Manuscripts/Fern Interference/Biotropica/Prototype Figures/Dry Ratio x BG.png", width = 70, height = 80, units = 'mm', res = 300)
GG_Drati_BG
dev.off()
## quartz_off_screen
## 2
tgc <- summarySE(DryMassExp, measurevar="LeafAveMass",
groupvars=c("Species"), na.rm = TRUE)
levels(tgc$Species) <- list("SASA" = "Sapindus saponaria",
"CYBA" = "Cymbopetalum baillonii")
GG_DLM_Spp <- ggplot(tgc, aes(x=Species, y=LeafAveMass, fill=Species)) +
geom_bar(stat = "identity", color = "black", alpha = .5) +
geom_errorbar(aes(ymin=LeafAveMass-se, ymax=LeafAveMass+se), width=.1) +
geom_text(label = c("b·", "a·"), aes(y = LeafAveMass+se, x = Species),vjust = -0.5, size = 5, family = "serif") +
scale_fill_manual(values=c("darkslategray3", "coral2")) +
ylab("\nDry Leaf Mass (g)") +
ylim(0, 0.18) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
text = element_text(family = "serif"),
axis.title = element_text(size = rel(1)),
axis.text = element_text(size = rel(0.9)), axis.ticks = element_line(color = "black", size = 0.8),
panel.border = element_rect(color = "black", fill = NA, size = 0.8),
strip.background = element_rect(color = "black", fill = "white", size = 0.8),
panel.grid = element_blank())
GG_DLM_Spp
png("~/Desktop/Universidad Javeriana/Manuscripts/Fern Interference/Biotropica/Prototype Figures/Dry Leaf Mass x Species.png", width = 70, height = 80, units = 'mm', res = 300)
GG_DLM_Spp
dev.off()
## quartz_off_screen
## 2
tgc
## Species N LeafAveMass sd se ci
## 1 CYBA 19 0.05157895 0.06784485 0.01556468 0.03270017
## 2 SASA 31 0.13354839 0.16817745 0.03020556 0.06168799
(0.0515789 - 0.1335484) / 0.0515789 * 100
## [1] -158.9206
tgc <- summarySE(DryMassExp, measurevar="LeafAveMass",
groupvars=c("Aboveground", "Species"), na.rm = TRUE)
tgc$Species_Aboveground <- paste(tgc$Species, tgc$Aboveground, sep = "_")
levels(tgc$Aboveground) <- list("AG C." = "Covered",
"No AG C." = "Exposed")
tgc$Species <- relevel(tgc$Species, ref = "Sapindus saponaria")
levels(tgc$Species) <- list("Sapindus" = "Sapindus saponaria",
"Cymbopetalum" = "Cymbopetalum baillonii")
tgc$Aboveground <- relevel(tgc$Aboveground, ref = "AG C.")
GG_DLM_AG <- ggplot(tgc, aes(x=Aboveground, y=LeafAveMass, fill=Species_Aboveground)) +
geom_bar(stat = "identity", color = "black", alpha = .5) +
geom_errorbar(aes(ymin=LeafAveMass-se, ymax=LeafAveMass+se), width=.1) +
geom_text(label = c(" b", "b", " a", "a"), aes(y = LeafAveMass+se, x = Aboveground),vjust = -0.5, size = 5, family = "serif") +
scale_fill_manual(values = list("paleturquoise", "paleturquoise4","lightsalmon", "lightsalmon4")) +
ylab("\nDry Leaf Mass (g)") +
ylim(0, 0.3) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
text = element_text(family = "serif"),
axis.title = element_text(size = rel(1)),
axis.text = element_text(size = rel(0.9)), axis.ticks = element_line(color = "black", size = 0.8),
panel.border = element_rect(color = "black", fill = NA, size = 0.8),
strip.background = element_rect(color = "black", fill = "white", size = 0.8),
panel.grid = element_blank()) +
facet_wrap(~Species, scales = "free", strip.position = "top")
GG_DLM_AG
png("~/Desktop/Universidad Javeriana/Manuscripts/Fern Interference/Biotropica/Prototype Figures/Dry Leaf Mass x AG.png", width = 70, height = 80, units = 'mm', res = 300)
GG_DLM_AG
dev.off()
## quartz_off_screen
## 2
tgc
## Aboveground Species N LeafAveMass sd se ci
## 1 AG C. Cymbopetalum 13 0.03153846 0.01675617 0.004647325 0.01012565
## 2 AG C. Sapindus 15 0.03600000 0.02323790 0.006000000 0.01286872
## 3 No AG C. Cymbopetalum 6 0.09500000 0.11220517 0.045807569 0.11775210
## 4 No AG C. Sapindus 16 0.22500000 0.19397594 0.048493986 0.10336248
## Species_Aboveground
## 1 Cymbopetalum baillonii_Covered
## 2 Sapindus saponaria_Covered
## 3 Cymbopetalum baillonii_Exposed
## 4 Sapindus saponaria_Exposed
(0.0360000 - 0.2250000) / 0.0360000 * 100
## [1] -525
(0.0315385 - 0.0950000) / 0.0315385 * 100
## [1] -201.2191
tgc <- summarySE(DryMassExp, measurevar="LeafAveMass",
groupvars=c("Belowground", "Species"), na.rm = TRUE)
tgc$Species_Belowground <- paste(tgc$Species, tgc$Belowground, sep = "_")
levels(tgc$Belowground) <- list("BG C." = "Unprotected",
"No BG C." = "Protected")
tgc$Species <- relevel(tgc$Species, ref = "Sapindus saponaria")
levels(tgc$Species) <- list("Sapindus" = "Sapindus saponaria",
"Cymbopetalum" = "Cymbopetalum baillonii")
tgc$Belowground <- relevel(tgc$Belowground, ref = "BG C.")
GG_DLM_BG <- ggplot(tgc, aes(x=Belowground, y=LeafAveMass, fill=Species_Belowground)) +
geom_bar(stat = "identity", color = "black", alpha = .5) +
geom_errorbar(aes(ymin=LeafAveMass-se, ymax=LeafAveMass+se), width=.1) +
geom_text(label = c("a", "a", "a", "a"), aes(y = LeafAveMass+se, x = Belowground),vjust = -0.5, size = 5, family = "serif") +
scale_fill_manual(values = list("skyblue4", "skyblue","mistyrose4", "mistyrose")) +
ylab("\nDry Leaf Mass (g)") +
ylim(0, 0.25) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
text = element_text(family = "serif"),
axis.title = element_text(size = rel(1)),
axis.text = element_text(size = rel(0.9)), axis.ticks = element_line(color = "black", size = 0.8),
panel.border = element_rect(color = "black", fill = NA, size = 0.8),
strip.background = element_rect(color = "black", fill = "white", size = 0.8),
panel.grid = element_blank()) +
facet_wrap(~Species, scales = "free", strip.position = "top")
GG_DLM_BG
png("~/Desktop/Universidad Javeriana/Manuscripts/Fern Interference/Biotropica/Prototype Figures/Dry Leaf Mass x BG.png", width = 70, height = 80, units = 'mm', res = 300)
GG_DLM_BG
dev.off()
## quartz_off_screen
## 2
tgc <- summarySE(DryMassExp, measurevar="Moisture",
groupvars=c("Species"), na.rm = TRUE)
levels(tgc$Species) <- list("SASA" = "Sapindus saponaria",
"CYBA" = "Cymbopetalum baillonii")
GG_Moist_Spp <- ggplot(tgc, aes(x=Species, y=Moisture, fill=Species)) +
geom_bar(stat = "identity", color = "black", alpha = .5) +
geom_errorbar(aes(ymin=Moisture-se, ymax=Moisture+se), width=.1) +
geom_text(label = c("b", "a"), aes(y = Moisture+se, x = Species),vjust = -0.5, size = 5, family = "serif") +
scale_fill_manual(values=c("darkslategray3", "coral2")) +
ylab("\nMoisture Content (%)") +
ylim(0, 60) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
text = element_text(family = "serif"),
axis.title = element_text(size = rel(1)),
axis.text = element_text(size = rel(0.9)), axis.ticks = element_line(color = "black", size = 0.8),
panel.border = element_rect(color = "black", fill = NA, size = 0.8),
strip.background = element_rect(color = "black", fill = "white", size = 0.8),
panel.grid = element_blank())
GG_Moist_Spp
png("~/Desktop/Universidad Javeriana/Manuscripts/Fern Interference/Biotropica/Prototype Figures/Moisture Content x Species.png", width = 70, height = 80, units = 'mm', res = 300)
GG_Moist_Spp
dev.off()
## quartz_off_screen
## 2
tgc <- summarySE(DryMassExp, measurevar="Moisture_Leaves",
groupvars=c("Species"), na.rm = TRUE)
levels(tgc$Species) <- list("SASA" = "Sapindus saponaria",
"CYBA" = "Cymbopetalum baillonii")
GG_MoistLeaves_Spp <- ggplot(tgc, aes(x=Species, y=Moisture_Leaves, fill=Species)) +
geom_bar(stat = "identity", color = "black", alpha = .5) +
geom_errorbar(aes(ymin=Moisture_Leaves-se, ymax=Moisture_Leaves+se), width=.1) +
geom_text(label = c("b", "a"), aes(y = Moisture_Leaves+se, x = Species),vjust = -0.5, size = 5, family = "serif") +
scale_fill_manual(values=c("darkslategray3", "coral2")) +
ylab("\nMoisture Content (%)") +
ylim(0, 60) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
text = element_text(family = "serif"),
axis.title = element_text(size = rel(1)),
axis.text = element_text(size = rel(0.9)), axis.ticks = element_line(color = "black", size = 0.8),
panel.border = element_rect(color = "black", fill = NA, size = 0.8),
strip.background = element_rect(color = "black", fill = "white", size = 0.8),
panel.grid = element_blank())
GG_MoistLeaves_Spp
tgc <- summarySE(DryMassExp, measurevar="Moisture",
groupvars=c("Aboveground", "Species"), na.rm = TRUE)
tgc$Species_Aboveground <- paste(tgc$Species, tgc$Aboveground, sep = "_")
levels(tgc$Aboveground) <- list("AG C." = "Covered",
"No AG C." = "Exposed")
tgc$Species <- relevel(tgc$Species, ref = "Sapindus saponaria")
levels(tgc$Species) <- list("Sapindus" = "Sapindus saponaria",
"Cymbopetalum" = "Cymbopetalum baillonii")
tgc$Aboveground <- relevel(tgc$Aboveground, ref = "AG C.")
GG_Moist_AG <- ggplot(tgc, aes(x=Aboveground, y=Moisture, fill=Species_Aboveground)) +
geom_bar(stat = "identity", color = "black", alpha = .5) +
geom_errorbar(aes(ymin=Moisture-se, ymax=Moisture+se), width=.1) +
geom_text(label = c("b", "b", "a", "a"), aes(y = Moisture+se, x = Aboveground),vjust = -0.5, size = 5, family = "serif") +
scale_fill_manual(values = list("paleturquoise", "paleturquoise4","lightsalmon", "lightsalmon4")) +
ylab("\nMoisture Content (%)") +
ylim(0, 65) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
text = element_text(family = "serif"),
axis.title = element_text(size = rel(1)),
axis.text = element_text(size = rel(0.9)), axis.ticks = element_line(color = "black", size = 0.8),
panel.border = element_rect(color = "black", fill = NA, size = 0.8),
strip.background = element_rect(color = "black", fill = "white", size = 0.8),
panel.grid = element_blank())
GG_Moist_AG
png("~/Desktop/Universidad Javeriana/Manuscripts/Fern Interference/Biotropica/Prototype Figures/Moisture Content x AG.png", width = 70, height = 80, units = 'mm', res = 300)
GG_Moist_AG
dev.off()
## quartz_off_screen
## 2
tgc <- summarySE(DryMassExp, measurevar="Moisture_Leaves",
groupvars=c("Aboveground", "Species"), na.rm = TRUE)
tgc$Species_Aboveground <- paste(tgc$Species, tgc$Aboveground, sep = "_")
levels(tgc$Aboveground) <- list("AG C." = "Covered",
"No AG C." = "Exposed")
tgc$Species <- relevel(tgc$Species, ref = "Sapindus saponaria")
levels(tgc$Species) <- list("Sapindus" = "Sapindus saponaria",
"Cymbopetalum" = "Cymbopetalum baillonii")
tgc$Aboveground <- relevel(tgc$Aboveground, ref = "AG C.")
GG_MoistLeaves_AG <- ggplot(tgc, aes(x=Aboveground, y=Moisture_Leaves, fill=Species_Aboveground)) +
geom_bar(stat = "identity", color = "black", alpha = .5) +
geom_errorbar(aes(ymin=Moisture_Leaves-se, ymax=Moisture_Leaves+se), width=.1) +
geom_text(label = c("b", "a", "a", "a"), aes(y = Moisture_Leaves+se, x = Aboveground),vjust = -0.5, size = 5, family = "serif") +
scale_fill_manual(values = list("paleturquoise", "paleturquoise4","lightsalmon", "lightsalmon4")) +
ylab("\nMoisture Content (%)") +
ylim(0, 65) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
text = element_text(family = "serif"),
axis.title = element_text(size = rel(1)),
axis.text = element_text(size = rel(0.9)), axis.ticks = element_line(color = "black", size = 0.8),
panel.border = element_rect(color = "black", fill = NA, size = 0.8),
strip.background = element_rect(color = "black", fill = "white", size = 0.8),
panel.grid = element_blank()) +
facet_wrap(~Species, scales = "free", strip.position = "top")
GG_MoistLeaves_AG
tgc <- summarySE(DryMassExp, measurevar="Moisture",
groupvars=c("Belowground", "Species"), na.rm = TRUE)
tgc$Species_Belowground <- paste(tgc$Species, tgc$Belowground, sep = "_")
levels(tgc$Belowground) <- list("BG C." = "Unprotected",
"No BG C." = "Protected")
tgc$Species <- relevel(tgc$Species, ref = "Sapindus saponaria")
levels(tgc$Species) <- list("Sapindus" = "Sapindus saponaria",
"Cymbopetalum" = "Cymbopetalum baillonii")
tgc$Belowground <- relevel(tgc$Belowground, ref = "BG C.")
GG_Moist_BG <- ggplot(tgc, aes(x=Belowground, y=Moisture, fill=Species_Belowground)) +
geom_bar(stat = "identity", color = "black", alpha = .5) +
geom_errorbar(aes(ymin=Moisture-se, ymax=Moisture+se), width=.1) +
geom_text(label = c(" b·", "a", " a·", "a"), aes(y = Moisture+se, x = Belowground),vjust = -0.5, size = 5, family = "serif") +
scale_fill_manual(values = list("skyblue4", "skyblue","mistyrose4", "mistyrose")) +
ylab("\nMoisture Content (%)") +
ylim(0, 65) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
text = element_text(family = "serif"),
axis.title = element_text(size = rel(1)),
axis.text = element_text(size = rel(0.9)), axis.ticks = element_line(color = "black", size = 0.8),
panel.border = element_rect(color = "black", fill = NA, size = 0.8),
strip.background = element_rect(color = "black", fill = "white", size = 0.8),
panel.grid = element_blank()) +
facet_wrap(~Species, scales = "free", strip.position = "top")
GG_Moist_BG
png("~/Desktop/Universidad Javeriana/Manuscripts/Fern Interference/Biotropica/Prototype Figures/Moisture Content x BG.png", width = 70, height = 80, units = 'mm', res = 300)
GG_Moist_BG
dev.off()
## quartz_off_screen
## 2
tgc <- summarySE(DryMassExp, measurevar="Moisture_Leaves",
groupvars=c("Belowground", "Species"), na.rm = TRUE)
tgc$Species_Belowground <- paste(tgc$Species, tgc$Belowground, sep = "_")
levels(tgc$Belowground) <- list("BG C." = "Unprotected",
"No BG C." = "Protected")
tgc$Species <- relevel(tgc$Species, ref = "Sapindus saponaria")
levels(tgc$Species) <- list("Sapindus" = "Sapindus saponaria",
"Cymbopetalum" = "Cymbopetalum baillonii")
tgc$Belowground <- relevel(tgc$Belowground, ref = "BG C.")
GG_MoistLeaves_BG <- ggplot(tgc, aes(x=Belowground, y=Moisture_Leaves, fill=Species_Belowground)) +
geom_bar(stat = "identity", color = "black", alpha = .5) +
geom_errorbar(aes(ymin=Moisture_Leaves-se, ymax=Moisture_Leaves+se), width=.1) +
geom_text(label = c(" b·", "a", " a·", "a"),
aes(y = Moisture_Leaves+se, x = Belowground),vjust = -0.5, size = 5, family = "serif") +
scale_fill_manual(values = list("skyblue4", "skyblue", "mistyrose4", "mistyrose")) +
ylab("\nMoisture Content (%)") +
ylim(0, 65) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
text = element_text(family = "serif"),
axis.title = element_text(size = rel(1)),
axis.text = element_text(size = rel(0.9)), axis.ticks = element_line(color = "black", size = 0.8),
panel.border = element_rect(color = "black", fill = NA, size = 0.8),
strip.background = element_rect(color = "black", fill = "white", size = 0.8),
panel.grid = element_blank()) +
facet_wrap(~Species, scales = "free", strip.position = "top")
GG_MoistLeaves_BG
tgc <- summarySE(DryMassExp, measurevar="Moisture_Leaves", groupvars=c("Treatment", "Species"), na.rm = TRUE)
tgc <- droplevels(subset(tgc, tgc$Species == "Cymbopetalum baillonii"))
levels(tgc$Treatment) <- list("None" = "No Competition",
"Only BG" = "Belowground Competition",
"Only AG" = "Aboveground Competition",
"Full" = "Full Competition")
tgc$Treatment <- relevel(tgc$Treatment, ref = "None")
GG_Moist_Int <- ggplot(tgc, aes(x=Treatment, y=Moisture_Leaves, fill=Treatment)) +
geom_bar(stat = "identity", color = "black", alpha = .5) +
geom_errorbar(aes(ymin=Moisture_Leaves-se, ymax=Moisture_Leaves+se), width=.1) +
geom_text(label = c(" a·", "a", "a", " b·"), aes(y = Moisture_Leaves+se, x = Treatment),
vjust = -0.5, size = 5, family = "serif") +
scale_fill_manual(values = list("skyblue4", "skyblue", "mistyrose4", "mistyrose")) +
ylab("\nMoisture Content (%)") +
xlab("Competition Treatment Combinations") +
ylim(0, 65) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
text = element_text(family = "serif"),
axis.title = element_text(size = rel(1)),
axis.text = element_text(size = rel(0.9)),
axis.ticks = element_line(color = "black", size = 0.8),
panel.border = element_rect(color = "black", fill = NA, size = 0.8),
strip.background = element_rect(color = "black", fill = "white", size = 0.8),
strip.text = element_text(face = "italic"),
panel.grid = element_blank()) +
facet_wrap(~Species, scales = "free", strip.position = "top")
GG_Moist_Int
png("~/Desktop/Universidad Javeriana/Manuscripts/Fern Interference/Journal of Applied Ecology/Review 1/Figures/CYBA_Moisture Content.png", width = 100, height = 80, units = 'mm', res = 300)
GG_Moist_Int
dev.off()
## quartz_off_screen
## 2
png("~/Desktop/Universidad Javeriana/Manuscripts/Fern Interference/Journal of Applied Ecology/Review 1/Figures/CombinedSpecies.png", width = 112.5, height = 250, units = 'mm', res = 300)
ggarrange(
GG_Est_Spp + labs(x = NULL),
GG_Height_Spp + labs(x = NULL),
GG_Diam_Spp + labs(x = NULL),
GG_LeafCount_Spp + labs(x = NULL),
GG_DMass_Spp + labs(x = NULL),
GG_DLM_Spp + labs(x = NULL),
GG_Drati_Spp + labs(x = NULL),
GG_MoistLeaves_Spp + labs(x = NULL),
align='h',
common.legend = F, nrow = 4, ncol = 2,
labels = c("a", "b", "c", "d", "e", "f", "g", "h")
)
dev.off()
## quartz_off_screen
## 2
png("~/Desktop/Universidad Javeriana/Manuscripts/Fern Interference/Journal of Applied Ecology/Review 1/Figures/Figure 2.png", width = 215, height = 250, units = 'mm', res = 300)
ggarrange(
GG_Est_AG + labs(x = NULL) + theme(strip.text = element_text(face = "italic")),
GG_Height_AG + labs(x = NULL)+ theme(strip.text = element_text(face = "italic")),
GG_Diam_AG + labs(x = NULL)+ theme(strip.text = element_text(face = "italic")),
GG_LeafCount_AG + labs(x = NULL)+ theme(strip.text = element_text(face = "italic")),
GG_DMass_AG + labs(x = NULL)+ theme(strip.text = element_text(face = "italic")),
GG_DLM_AG + labs(x = NULL)+ theme(strip.text = element_text(face = "italic")),
GG_Drati_AG + labs(x = NULL)+ theme(strip.text = element_text(face = "italic")),
GG_MoistLeaves_AG + labs(x = NULL)+ theme(strip.text = element_text(face = "italic")),
align='h',
common.legend = F, nrow = 4, ncol = 2,
labels = c("a", "b", "c", "d", "e", "f", "g", "h")
)
dev.off()
## quartz_off_screen
## 2
png("~/Desktop/Universidad Javeriana/Manuscripts/Fern Interference/Journal of Applied Ecology/Review 1/Figures/Figure 3.png", width = 215, height = 250, units = 'mm', res = 300)
ggarrange(
GG_Est_BG + labs(x = NULL)+ theme(strip.text = element_text(face = "italic")),
GG_Height_BG + labs(x = NULL)+ theme(strip.text = element_text(face = "italic")),
GG_Diam_BG + labs(x = NULL)+ theme(strip.text = element_text(face = "italic")),
GG_LeafCount_BG + labs(x = NULL)+ theme(strip.text = element_text(face = "italic")),
GG_DMass_BG + labs(x = NULL)+ theme(strip.text = element_text(face = "italic")),
GG_DLM_BG + labs(x = NULL)+ theme(strip.text = element_text(face = "italic")),
GG_Drati_BG + labs(x = NULL)+ theme(strip.text = element_text(face = "italic")),
GG_MoistLeaves_BG + labs(x = NULL)+ theme(strip.text = element_text(face = "italic")),
align='h',
common.legend = F, nrow = 4, ncol = 2,
labels = c("a", "b", "c", "d", "e", "f", "g", "h")
)
dev.off()
## quartz_off_screen
## 2